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Arcjet thrusters are being actively considered for use in Earth orbit maneuvering 
applications. Satellite station-keeping is an example of a maneuvering application requiring 
the low thrust, high specific impulse of an arcjet. Experimental studies are currently the 
chief means of determining an optimal thruster configuration. Earlier numerical studies 
have failed to include all of the effects found in typical arcjets including complex 
geometries, viscosity and swirling flow. 

Arcjet geometries are large area ratio converging-diverging nozzles with 
centerbodies in the subsonic portion of the nozzle. The nozzle walls serve as the anode 
while the centerbody functions as the cathode. Viscous effects are important because the 
Reynolds number, based on the throat radius, is typically less than 1,000. Experimental 
studies have shown a swirl or circumferential velocity component stabilizes a constricted 
arc. 

This dissertation describes the equations governing flow through a constricted arcjet 
thruster. An assumption the flowfield is in local thermodynamic equilibrium leads to a 
single fluid plasma temperature model. An order of magnitude analysis reveals the 
governing fluid mechanics equations are uncoupled from the electromagnetic field 
equations. A numerical method is developed to solve the governing fluid mechanics 
equations, the Thin Layer Navier-Stokes equations. A coordinate transformation is 
employed in deriving the governing equations to simplify the application of boundary 
conditions in complex geometries. 

An axisymmetric formulation is employed to include the swirl velocity component 
as well as the axial and radial velocity components. The numerical method is an implicit 
finite-volume technique and allows for large time steps to reach a converged steady-state 
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solution. The inviscid fluxes are flux-split and Gauss-Seidel line relaxation is used to 
accelerate convergence. 

Converging-diverging nozzles with exit-to-throat area ratios up to 100:1 and 
annular nozzles were examined. Comparisons with experimental data and previous 
numerical results were in excellent agreement. Quantities examined included Mach number 
and static wall pressure distributions, and oblique shock structures. As the level of swirl 
and viscosity in the flowfield increased the mass flow rate and thrust decreased. The 
technique was used to predict the flow through a typical arcjet thruster geometry. Results 
indicate swirl and viscosity play an important role in the complex geometry of an arcjet by 
shifting the Mach contours upstream and reducing the mass flow rate and thrust. 
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Chapter 1 


Introduction 


1.1 Background 

An arcjet is an electrothermal propulsion device in which a gas propellant is heated by 
means of a high temperature arc. The propellant expands through a large area ratio nozzle to 
produce thrust. A typical constricted arcjet configuration is shown in Fig. 1.1. The 
propellant gas enters a chamber, flows over and past an electrically conducting rod serving 
as a cathode. An electric arc is established between the cathode and the nozzle walls, which 
serve as an anode. The arc heats the gas to a high temperature, causing an expansion and 
resultant acceleration of the gas out through a constant area constrictor region and into the 
divergent portion of the nozzle where the arc attaches diffusely. 



Figure 1.1 Constricted Arcjet Schematic 

The constrictor region stabilizes the inherently unstable columnar arc. A "cool" 
boundary layer surrounds the high temperature (» 10,000 K) arc, preventing erosion of the 
nozzle walls. The tangential or swirl velocity component imparted to the inlet gas flow also 
provides a stabilizing force by establishing a radial pressure gradient. The pressure gradient 
opposes any kinking occurring in the arc column [1]. 

The arcjet was identified as a prime candidate for space missions requiring high 
specific impulse (1,000 -2,000 seconds) and low thrust (* llbf) propulsion devices in the 
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1960’s [2]. In contrast the maximum specific impulse of chemical rockets is on the order of 
200-400 seconds, and is limited by the chemistry of available fuels [3,4]. A typical mission 
requiring a high specific impulse would be to maintain a satellite's altitude in a geo- 
stationary orbit. Arcjets requiring 1-200 kilowatts (Kw) of power were being developed for 
a mission such as this in the 1960's. The research and development effort had led to 
continuous operation of up to 500 hours and overall efficiencies of up to 55 percent. 
However, excessive weight in the electric power generating systems was cited as a main 
development hurdle. 

Experimental development efforts have identified several items as greatly affecting 
arcjet performance. As mentioned previously, imparting a swirl or tangential velocity 
component to the inlet gas flow has been found (as early as 1909, [5]) to stabilize the arc 
column and reduce erosion of the cathode and nozzle walls. Recent work at NASA by 
Curran [6] has confirmed the need for a strongly swirling flow. 

Unsuprisingly, the cathode shape, nozzle geometry, and constrictor length were 
also found [6] to have an affect on gross performance parameters such as specific impulse 
and thrust as well as on localized phenomenon such as the downstream current attachment 
location, nozzle wall temperatures, etc. In addition, the operating conditions and 
dimensions of arcjet thrusters lead to Reynolds numbers below 1,000, indicating the flows 
are laminar yet highly viscous [7]. Any constricted arcjet model developed needs to include 
the effects of swirl, geometry, and viscosity. 

1.2 Previous Research 

Because of the modest computational facilities available in the 1960's, the research 
effort regarding arcjets was largely experimental. The earliest modeling efforts were 
applicable in a constant area region downstream of the cathode tip (See Fig. 1.2). Among 
the earliest modeling efforts were those of Stine and Watson [8] and Watson and Pegot [9]. 
Their simplifying assumptions of unswirled, fully developed, one dimensional flow led to 
an uncoupled energy equation which accurately modeled the enthalpy distribution in the 
arc. 

Later research by Neuberger [10,11] and Shaeffer [1] involved an axisymmetric 
model in a constant area geometry, with viscous, swirling flow. Shaeffer also included 
turbulence and radiation transport effects in his model. Neuberger's results indicate the 
effect of swirl was small in geometries with large constrictor length/diameter ratios but in 
small ratio geometries higher averaged enthalpy and efficiencies were reached. Shaeffer 
demonstrated the stabilizing effect of swirl on the arc column decreased unless secondary 
gas injection occurs. 

However, in examining constant area geometries both Neuberger and Shaeffer 
assumed radial velocities were small compared to the axial and tangential (swirl) velocities 
and axial gradients were negligible in comparison to radial gradients. These assumptions 
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are essentially boundary layer assumptions (also called the quasi-cylindrical assumption) 
and when applied to an axisymmetric, swirling flow are termed the weak swirl 
approximation [13]. In this context, they result in a greatly reduced radial momentum 
equation 


pw 2 _ 3P 


which is a balance between the centrifugal force and the radial pressure gradient. When no 
assumptions are made about radial velocity or axial gradient magnitudes a strongly swirling 
flow can be modeled and consists of the axisymmetric Navier-Stokes equations. The radial 
momentum equation for a strongly swirling flow is given by 


1 a i ,\ . a , „ \ P w 2 5P . i 3 \ . » /, \ . '■ee 

r 5r( r P v ) + &( rpuv ) - — = + r dF ( + & (t " ) + ~ 


x ee 


Obviously a great deal of the physics of the flow is neglected in an assumption of weakly 
swirling flow. Yet in a constricted arcjet, in which the geometry varies greatly with axial 
location, no assumptions can be made about radial velocity or axial gradient magnitudes. 
Thus, an accurate treatment of the flow from a fluid mechanics standpoint requires a 
strongly swirling flow to be modeled using the Navier-Stokes equations. 


Tangential Row 


Cathode 



Figure 1.2 Characteristic regions of flow through a wall and 
swirl stabilized arcjet thruster (after Neuberger, 
1975) 
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It should also be mentioned that a strongly swirling flow is generally indicated by 
an inlet swirl number, Sj , greater than 1 .0. Sj is defined by 


Si 



puwr 2 dr 



pu 2 r dr 


which is seen to be the axial flux of angular momentum divided by the inlet wall radius 
times the axial flux of axial momentum, and is a direct measure of the level of swirl at the 
nozzle inlet. The Neuberger and Shaeffer studies were limited to swirl numbers less than 
1.0 and thus were limited to examining the effects of weakly swirling flow. In general, 
weakly swirling and unswirled flows as well as strongly swirled flows can be examined 
using the axisymmetric Navier-Stokes equations. 

More recently, Nishida, et. al [13] , modeled a realistic, axisymmetric arcjet 
geometry, i.e., a diverging nozzle exit as in Fig. 1.1, assuming unswirled, inviscid flow. 
Also included were nonequilibrium ionization-recombination effects for an argon 
propellant. The mass flow rate, discharge current, and geometry were that of an 
experimental thruster undergoing testing. Numerical results indicated a decrease in thrust 
and specific impulse as the constrictor length increased. No comparisons with experimental 
data were provided. 

In related work, [14,15], researchers at the University of Stuttgart have been 
developing numerical models of magnetoplasmadynamic (MPD) thrusters (in which high 
current densities generate self-magnetic fields producing substantial electromagnetic 
acceleration of the flow and hence greatly complicate the governing equations). Their initial 
models have been inviscid and partly two-dimensional in nature. Predicted current density 
distributions in cylindrical geometries have compared favorably with experimental results. 
Further development efforts are proceeding to reduce the mass flow errors (up to 10%) and 
implement more realistic subsonic inlet boundary conditions in place of the supersonic inlet 
conditions currently used. 

Worthy of special mention are the detailed analyses of the anode contraction region 
in high intensity arcs operating at atmospheric pressure performed by Pfender and his 
students, [16-19]. The analyses used a two temperature model, which accounts for 
nonequilibrium effects near the anode or cathode surface by including energy equations for 
electrons and heavy particles (ions and neutral atoms) in addition to the continuity and 
momentum equations. 

No attempt has been made in any of the arcjet modeling efforts discussed above to 
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include a detailed analysis of regions near the cathode tip or anode wall. The complexity of 
an analysis of the entire flowfield precludes the ability to accurately predict details in these 
regions of tremendous current, temperature, and velocity gradients. Eventually 
electrothermal thruster models would need to incorporate accurate boundary conditions in 
these regions. 

In summary, none of the works examined to date have incorporated the ability to 
model the effects of strongly swirling, viscous flows throughout an entire arcjet thruster. 

1.3 Present Research Effort 

In this dissertation a numerical model of a constricted arcjet thruster is developed. 
The formulation includes all of the flow features found to have significant effects on the 
operation of a thruster including realistic geometries, viscous effects and a strongly 
swirling flow. 

The continuity, momentum, energy, and electromagnetic field equations are derived 
in general terms assuming the flowfield is in local thermodynamic equilibrium (LTE). Thus 
the flow is assumed to have a single mean temperature (electrons, ions and neutral atoms 
have the same temperature) at any point. This assumption is valid everywhere except in 
regions close to the anode and cathode surfaces. 

An order of magnitude analysis demonstrates the fluid mechanics equations 
(continuity, momentum, and energy) are uncoupled from the electromagnetic field 
equations when insignificant terms are neglected. The only link between the two equation 
sets is through the electrical conductivity, which is a highly non-linear function of pressure 
and temperature. 

Both equation sets are formulated in transformed coordinates. The fluid mechanics 
equations solved are the Thin Layer Navier-Stokes (TLNS) equations. MacCormack's [20] 
implicit Gauss-Seidel Line Relaxation method is modified to solve the axisymmetric form 
of the TLNS equations. The iterative technique needed to obtain a converged solution of the 
TLNS algorithm and the electromagnetic field equations (solved by a forward in time, 
centered in space algorithm) is also discussed. 

A grid generation algorithm is also developed to create the numerical mesh. The 
algorithm allows grid to be clustered near upper and lower surfaces as well as at the throat. 
Thus, regions containing large gradients such as in viscous boundary layers or in the 
constrictor can be adequately resolved. In addition, the algorithm can generate mesh 
orthogonal to upper and lower surfaces, allowing the effect of grid orthogonality on 
solution accuracy to be determined. 

The TLNS portion of the formulation (valid for flow through a thruster with the arc 
turned off, i.e., cold flow) is validated by comparison with experimental data obtained for 
flows through converging-diverging nozzles. Finally, the effect of swirl on flow through a 
geometry representative of a thruster is examined. 
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Chapter 2 


Governing Equation Formulation 


2.1 Introduction 

In this chapter the equations governing the flow through a constricted arcjet thruster are 
presented and non-dimensionalized. An order of magnitude analysis is performed to 
determine which terms in the governing equations have a negligible effect on the thruster 
flowfield. A non-dimensionalization is applied to the reduced set of governing equations. 
Finally, a coordinate transformation is applied to derive the governing fluid mechanics and 
electromagnetic field equations in transformed coordinates. 

2.2 Basic Assumptions 

The following assumptions will be made in order to derive the governing equations 
for a vortex stabilized arcjet thruster: 

1) The flowfield is assumed to be in local thermodynamic equilibrium (LTE). Thus 
the electron, ion, and heavy or neutral particle temperatures are assumed equal and can be 
equated to a single mean plasma temperature. This results in a single fluid plasma model. 
This continuum approach is not valid near the anode and cathode surfaces where kinetic 
theory, i.e. a two temperature plasma model, should be used, [1, 10, 16]. 

2) The plasma is assumed optically thin, i.e., reabsorption of emitted radiation is 
neglected, [10, 19]. 

3) The net space charge density, p e , is assumed to be zero. This results in a quasi- 
neutral plasma. Appendix A examines the validity of this assumption, [10, 17]. 

4) External magnetic fields and the effects of gravity arc neglected. 

5) The Lorentz force (self induced magnetic force) is included in the momentum 
equation and the transport of electron enthalpy due to the drift (diffusion) of the electrons is 
included in the energy equation. These terms are, respectively [16, 18] 
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jxB ; ^fcj VT 
2 e 

Expressions for energy change due to volume expansion and viscous dissipation 
are also included in the energy equation. An order of magnitude analysis will be performed 
to determine what terms are retained in the governing equations. 

Using the above five assumptions the governing equations for the plasma may be 
written [21, 22, 23]. These include the continuity equation 


the momentum equation 


^ + V.pv = 0 
dt 


p Dy. = - VP - V-x + j x B 
K Dt 


( 2 . 1 ) 

( 2 . 2 ) 


where t is the shear stress tensor and if Stoke's assumption for the bulk viscosity is made, 
may be expressed as 


x 


ij 



fixj dxi) 3 ° Xk 


The conservation of energy equation takes the following form 


*2 

pDej. = _ v * . ypv . V-(x • v) + J - - S R + ^ j V T 

r Dt M v ' c 2 e 

or (2.3) 

^ + V E, v = -V q - V Pv - V-(t v) + £• - Sr + j V T 
dt 

where 

q = - kVT 

e t = total energy per unit mass 

Ej = pe, = total energy per unit volume 

kg = Boltzmann's constant 

e = unit of electric charge 

Cp = specific heat at constant pressure/unit mass 
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The terms on the right hand side of the Eq. (2.3) represent a change in energy due 
to heat conduction, volume changes, viscous dissipation. Joule heating, radiation, and the 
electron enthalpy flux. 

Neglecting electron pressure effects and ion slip [22, 23] the Generalized Ohm’s 
Law can be written 


♦ 

^ — +j = c(e + vxb) - p^jxB (2.4) 

He e 2 a t Dee 

where 

a = electrical conductivity 
% = electron number density 
1 % = electron mass 
e = unit of electric charge 

The three terms on the right hand side of the last equation represent the applied 
electric field, the induced field due to the bulk movement of the plasma, and the Hall effect 
due to the interaction of the current with the magnetic field, respectively. 

In addition the electromagnetic fields obey Maxwell's equations, which are, for a 
net space charge density equal to zero ( Assumption 2); the induction law 

VxE = -— (2.5) 

at 


This implies V • B = 0 as may be seen by taking the divergence of Eq. (2.5) (since the 
divergence of the curl of a vector field is zero). 

By neglecting the magnetization of the plasma, the magnetic permeability is 
essentially that of a vacuum, and the current law is 


where 


V x B = n Mo j (2.6) 

|i Mo = permeability of vacuum 
= 4 7t 10-7 henry/m 


If the divergence of the last equation is taken we obtain a current continuity 
equation: 

V • j = 0 


8 


In addition we will assume the plasma obeys the perfect gas law (even though 
ionized gases should not be considered perfect) 

P = p R T (2.7) 


where R = R / M.W. 

R = universal gas constant 
M.W. = molecular weight 

Finally, we have the following set of dimensional governing equations (where the A 
indicates a dimensional quantity): 


Continuity Eq 


dp C. _ 
- 7 *+ V.pv = 0 

dt 


Ar 

Momentum Eq 


AN A * A 

A A A ♦ ^ 

p Dy. = . VP - V T + j x B 
Dt 


Energy Eq 


dE t 

dt 


+ 


As 

V -E t v 



- V q - V -Pv - 




Sr + 


ika/f 

2 a J 

e 


VT 


Generalized Ohm's Law 

Hr + j = of E + vxb| - ^rjxB 
n;e 2 dt \ / n * e 


Induction 


V 



dB 

A 

dt 


Current 


Equation of State 


VxB = p M oj 


V • B = 0 

A A 

V . j = 0 


P = p R T 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 
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2.3 Order of Magnitude Analysis 


The total number of equations including the equation of state is 15. The 22 
dimensional variables are (where the superscript " A " indicates a dimensional value of 
the quantity) 


v, B, E, j, p, P, T, o, k, p, n t , c^, Sr, y 


Letting: (where the subscript " 0 " indicates a dimensional reference value of 

the quantity) 


V = V 

• Vo 

B = B • B 0 

E = E • E 0 

A 

♦ ♦ 
j=j • 

a 0 E 0 

P = P Po 

W 

tl 

o 

H> 

II 

H 

T 0 

a = o • o o 

k = k • k 0 

P = P 

' Po 

n e = n e • n eo 

Cp = Cp ■ Cpo 

Sr = : 

Sr • Sr 0 

Xi = Xi r th 

t = t (r* / v 0 ) 

A 

X = X 

• (Po • v 0 / rth) 


<1 > 
II 

< 

& 

D/Dt = 

= D/Dt ■ v 0 / rth 

E t = E t 

(Po ‘ Cpo * T 0 ) 


and where 


P 0 = Po R T 0 

R = C po - C vo 


To ~ c po I c vo 


and defining j 0 > E 0 , and B 0 with: 


jo — ^ E 0 — j c / O q 


- Pmo °° Eo Tq, 


We can now obtain the familiar non-dimensional parameters: Re, M, Pr, etc.. 


Re = Reynolds number = v 0 r^ p^ = Inertia forces / viscous forces 
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Mo = Mach number = v 0 /(y 0 R T 0 )l/2 = Fluid velocity / speed of sound 

Pr = Prandtl number = Cj*, |i 0 / k Q = Momentum diffusivity / thermal diffusivity 

Pm = Magnetic pressure number = B 0 ^/ ( P 0 )i Mo ) 

= Pressure of magnetic field forces / static pressure 

R l = Electric heating number = Go Eo^ r^ / ( ko-T 0 ) 

= Electric energy supplied / energy transport by heat conduction 

R s = radiation number = S^r^ / ( k^To ) 

= Radiated energy / energy transport by heat conduction 


Rc = conduction number = ko Oo Eo / ( e(ko / r^)) 

= electrical conduction / thermal conduction 

Rp = magnitude of current fluctuation = nvao-vo / n^-e^-r^ 

R h = Hall parameter = o 0 Bo / n^ e = Hall current / current parallel to electric field 


Rjn = Magnetic Reynolds number = v 0 -Bo/E 0 = Vlij^ OoTti, 

= induced electric field / applied electric field 

Substituting the above quantities into Eqs. (2.8) - (2.14) yields the following non- 
dimensional equations: 


Continuity Eq 


+ Vpv = 0 
9t 


(2.15) 


Momentum Eq p ^ — VP - V t + — j X B 

Yo Yo M^ 


(2.16) 


Energy Eq 


9 E t 

at 


+ V-E t "v = 


V-(kVT) - V-Pv - V-(t-v) 


Re Pr Yo Re 

+ . -BS-Sr + J -v T 

Re Pr a Re Pr R 2 Re Pr J 


(2.17) 
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Generalized Ohm's Law 



+ j = 0 E 


+ RjnOvxB - Rh n^j x B 


( 2 . 18 ) 


Induction Law V xE — — Rin (2.19) 

dt 

Current Law V X B = j (2.20) 

Equation of State P = pT (2.21) 

Thus, the dimensional values of: P 0 , p 0 , T 0 , \i D , c^, k 0 , S ro , a 0 , E 0 , r^,, y 0 , n eo 
are needed to estimate the similarity parameters. A review of recent arcjet literature [6, 24, 
25, 26] provides the following representative values for a pure nitrogen propellant gas in an 
arcjet: 

P 0 = 1 atm = 1 xlO 5 T 0 = 10,000 K 

r^ = 2.0 xlO' 4 m n^, = 2 xlO 21 electrons / m 3 

j 0 = 10 amp / ( 7t rjji 2 ) = 7.96 xlO 2 amp / m2 

Mass flow rates = 2.7 - 6.1 xlO'^ kg / sec 

The following values were extracted from the nitrogen equilibrium transport 
properties compiled by Liu [16] from 12 independent sources 

Cpo = 5.55 xlO 3 J / kg K p c = 1.68 xlO* 2 kg / m 3 

H 0 = 2.57xlO- 4 Ns/m 2 k 0 =1.18W/mK 

C 0 = 2.84 xlO 3 1 / ohm m S ro = 1.47 xlO 8 W / m 3 

Thus, the following values can now be calculated 

E 0 = jo / °o = 28,028 volt / m 
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Bo = M-m 0 o 0 E 0 % = 2.00 xlO' 2 kg / coul sec 
In addition, for nitrogen: 


R = HI M.W..= (8.314 J / mole K) / (28 xlO' 3 kg / mole) 
= 296.93 J/ kg K 
Yo = c po / (Cpo ' R) = 1-056 
c Q = (Yo R T 0 ) 1/2 = 1 J70 m / sec 


Values for the non-dimensional parameters independent of the velocity, v 0 i.e., Pr, 
P M , Rl, Rs> Rc> and Rh* can now be calculated. The remaining non-dimensional 
parameters : Mq, Re, R F , and R m , will be calculated for a range of velocities. Tables 2.1- 
2.4 contain the results of these calculations for nitrogen under the above conditions, with 
r^ = 2.0 xl0~4,and 2.0 xlO" 3 m. 

In order to neglect the effect of a term in a governing equation the magnitude of the 
non-dimensional parameters of the term need to be 1-2 orders less than the next largest term 
in the governing equation. Examination of the parameters within the Momentum equation, 
Table 2.1, shows the only term that might be neglected is the Lorentz force but it is not at 
least an order of magnitude less than the pressure and viscous terms. Of course, the effect 
of the self induced magnetic field on the flow can only be examined by retaining this term. 

Examination of the Energy equation. Table 2.2 shows the viscous dissipation term 
can be neglected. The radiation term could be neglected but will be retained. This term will 
appear as a source term in the governing equations when they are placed into the proper 
format used in the numerical procedure which will be discussed later. 

Examination of the parameters in Ohm's Law, Table 2.4, indicates that the unsteady 
current term and the induced electric field can be neglected. However, the Hall effect, 

which contains ajxB term, is almost of sufficient magnitude to retain the term along with 
the applied electric field, E, which has an order of magnitude of one. Yet, in the 
Momentum Equation the Lorentz force, which also contains j x B. is of a much lower 
magnitude. 

Examination of the non-dimensionalization reveals 


Lorentz Force => 


_Zm_ 

YoM^ 




Hall Effect => R H - ^ 


Induced field => Rin * f* 
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As the length parameter, r^, of the non-dimensionalized equations increases both 

the Lorentz force and Hall effect terms decrease as shown above. Table 2.3 reveals that if 

r^ is increased by a factor of ten the Hall effect and induced field terms are negligible when 

compared to the applied field. The induced field increases proportional to but it is of such 

a small magnitude that it may still be neglected. Tables 2.1 and 2.4 show that when r^ is 

* — ♦ 

increased by a factor of ten the j x B terms can both be neglected in their respective 
governing equations. 

Thus, the order of magnitude analysis demonstrates that for completeness when 
modeling an arcjet the Hall effect term would be retained in Ohm's Law if the Lorentz force 
was included in the Momentum Equation. In this study, as in other arcjet studies [13, 16, 
17] the effect of the induced field and Hall effect terms will be considered negligible and the 
Lorentz force will be retained in the Momentum Equation. 

The final expression to examine is Eq. (2.19), the Induction Law. R IN , whose 
magnitude is shown in Table 2.4, is sufficiently small to neglect the unsteady term. 

The resulting set of non-dimensional governing equations is 


Continuity Eq 


^ + V.pv = 0 
9t 


( 2 . 22 ) 


Momentum Eq p 1 — VP - V x + ^ j x B (2.23) 

01 Re YoK; 


Energy Eq 


1 + VE, v = 
dt 

— L_ V- kVT - ,0 V Pv 

Re Pr Yo 


+ ii 2 . 

J*S_ s R + 5-Jk- j V T 

Re Pr R 2 Re Pr J 

(2.24) 

RePr o 

Generalized Ohm's Law 

j = oE 

(2.25) 

Induction Law 

o 

II 

t w 

X 

> 

(2.26) 

Current Law 

V x B = j 

(2.27) 

Equation of State 

P = p T 

(2.28) 


14 


Table 2.1 


Momentum Equation Order of Magnitude Analysis 



TWEEMBtM 




msm 



HMW.ILBafiTTOI 

17.7 


m*w 

HKSStSSl 


uatmaam 

177 




| 


1,770 


mm 

H!£33S9i 

ihx&esbih 

■BS&lSHi 

3,541 


46. 




r,h= 2.0 

x 10‘ 3 m 


mwmrm 


Re 

Pressure Force 

Viscous Force 

Lorentz Force 

17.7 




■SQlSISi 


177 

0.10 

«*» 


0.43 E-01 

■Kfisissm 

1,770 

1.00 

WtIM 


wMsssssmtm 


3,541 

2.00 

460 

0.24 E+00 

0.22 E-02 

0.18 E-05 


— 1 — => Pressure Force 
To MS 


-J- => Viscous Force 
Re 


— M =* Lorentz Force 

7o m2 


P 


Dv= 1 — VP 

** Yo Mo 


. _L 


Re 


V-x 


+ — M ] x B 

Y<>mJ 
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Table 2.2 


Energy Equation Order of Magnitude Analysis 
for r th = 2.0 x 10*4 m 


1 rth- 2.0 x 10~ 4 m 1 


Thermal Conduction 


Viscous Dissipation 

17.7 

“ 0.36 E+01 

12212 ! 


177 




1,770 


0.53 E-01 

0.24 E-02 

3341 


0.53 E-01 



Joulean Heating 

Radiation 

Electron Enthalpy 

17.7 

~ 0.27 E+02 ' 


f— 0.42 E+01 

177 

" 0.27 E+01 



1,770 




3,541 

0.14 E+00 


0.21 E-01 


— *— => Thermal Conduction 

RePr 


Yo- 1 
Yo 


=> Volume Changes 


Mi(Yo-l) 

Viscous Dissipation 

Re 


r l 

RePr 


=> Joulean Heating 


— => Radiation 
RePr 


RePr 


Electron Enthalpy 


9 E t 
9t 


+ 


VE, v = 


— I— V-(kVT) - VPv 

RePr I I y 0 


m5(yo-i' 

Re 


+ 


Rl j 2 
RePr o 


Rs 

RePr 


Sr + 


5 Rc 
2 RePr 



T 
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Table 2.3 


Energy Equation Order of Magnitude Analysis 
for r th = 2.0 x 10'3 m 


| rth m 2.0 x 10 3 m | 

MKHl 

Thermal Conduction 

Volume Changes 

Viscous Dissipation | 

17.7 




177 

1 



17770 


»ki 


3,541 

0.18 E-02 

0.53 E-01 



Joulean Heating 


Electron Enthalpy 

17.7 

r 0.27 E-01 

0.22 E-01 

0.42 E-01 

111 



" 0.42 E-02 " 

1,770 

1 

0.22 E-03 


£541 



0.21 E-03 | 


— * — => Thermal Conduction 
RePr 


Yo ~ 1 
Yo 


=> Volume Changes 


n^(yo-i) _ , . 

=> Viscous Dissipation 

Re 


— => Joulcan Heating 
RePr 6 


=> Radiation 

RePr 


=> Electron Enthalpy 

RePr 


dEt 

at 


+ 


VE t v = 


— 1— V-(kVT) - - — - VPv 
Re Pr V / Yo 


m£(yo-i) 

Re 


+ _Rk_i! . -Rs- 

RePr o Re Pr 


Sr + 


5 Rc 
2 RePr 


j-V 


T 
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Table 2.4 


Ohm's Law Order of Magnitude Analysis 


| rth- 2.0 x 10’ 4 m | 


Unsteady 



17.7 

0.45 E-05 


0.18 E+00 

177 




1,770 * 

0.45 E-03 



3,541 ] 



■■KSE23SBM 


| r//,= 2.0 x 10’ 3 m 

mmr^nu 

Unsteady 

Induced Field 


17.7 



0.18 E-01 

177 

0.45 E-05 



TTHO 



0.18 E-01 

3,541 


0.25 E-01 

0.18 E-01 4 


Rp => Unsteady term 


Rin => Induced field Rh => Hall Effect 


Rf £ 


<Ll 

at 


+ j = oE 
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2.4 Non-dimensional Governing Equations 

The preceding non-dimensional analysis was useful in demonstrating that the 
viscous dissipation term in the energy equation and the induced field and Hall effect terms 
in Ohm's Law can be neglected. In order to provide a consistent set of notation with which 
to compare the governing equations and results of this investigation to the work of other 
researchers, another group of non-dimensional parameters will be introduced. 

The final set of non-dimensional parameters will be based on the dimensional inlet 
stagnation values of density, p 0 , temperature, T 0 , speed of sound, c D , kinematic viscosity, 
Po, and the throat radius, r^. The following non-dimensional quantities are defined, where 

the subscript " 0 " and the superscript " A " are as noted at the beginning of Section 2.3, 


t> 

11 

<t> 

• c 0 

B = B • B 0 

E = E • E 0 

A 

j=j • 

Go Eq 

P = P • Po 

P = P • (po -eg) 

T = T 

•T 0 

a = o • g 0 

k = k • ko 

) 

II 

T: 

• V Lo 

n e = n e * n eo 

Cp = Cp * Cpo 

Sr = 

Sr • Sr 0 

Xi = Xj T lh 

Cy = Cy * Cpo 

A 

X = X 

• (Po • Co / r lh ) 

*t = t (fth / c 0 ) 

E t = E t • (po • eg) 

D/D t 

= D/Dt • c 0 / r,h 

<3 > 
II 

<3 

c t = e t • eg 


The non-dimensional speed of sound, c, is 

c2 „ 2. = jlAJL = i T = 

cl 7. • R • T. 1° P 
where 7o = ^ «*1 

Cvo Cv t ' v 
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The non-dimensional perfect gas relationship becomes 


P = -LpT 

Yo 


• . yR r 

with <?=— Cv =^r 

~ (£/m.w.) 

where R = = - — - = Cp - c v 

Cp 0 Cpo 

R = universal gas constant 
M.W. = Molecular weight 

The non-dimensional total energy per unit volume is 

E, = p e, = p (e; + u? ^ . ±w ? ) 

where u, v, and w are the axial, radial, and circumferential velocity components, 
respectively, and ej is the non-dimensional internal energy per unit mass given by 

. = £ = Si = i_jh = _i__e 

' 3 3 (*-')* < Y ‘ 1 < P 

In addition, the pressure P can be expressed as 

p = p (y- 1) (<=, • “ 2 + y2 ±y2 ) 

Only one non-dimensional parameter, the Reynolds number, must be redefined as 

Re = Reynolds number = c 0 r* p 0 / |io = Inertia forces / viscous forces 

with the re main ing non-dimensional parameters unchanged. 

With the current non-dimensionalization the non-dimensional governing equations 

become 


Continuity Eq 


Vpv = 0 


( 2 . 29 ) 
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Momentum Eq 


(2.30) 


pEv = . VP - -L V x + 5 m. j x B 
K Dt Re Yo 


Energy Eq 


0 E t 


_L 


a , +VE ‘ v RePr^-l) 


V-(kVT) - V Pv 


, Rl r 

Rs 

i 5 | V T 

(2 31) 

Re Pr (y 0 - 1) a 

Re Pr (Yo 

-1) Sr 2RePr(Yo-l) J 


Generalized Ohm's Law 






j = oE 

(2.32) 

Induction Law 


o 

II 

T W 
X 
> 

(2.33) 

Current Law 


V x B = j 

(2.34) 

Equation of State 


P = -LpT 
Yo 

(2.35) 


Making the assumption the flow is axisymmetric ( — = 0 ) and representing the axial, 

90 

radial, and circumferential or swirl velocity components as u, v, and w, respectively, we can 
write Eqs. (2.29-2.31) in vector form as: 


0Q 0F' 0G’ _ 

— + — + + H= 0 

0t 0z 0r 


(2.36). 


where the rows in the vectors Q, F, G’, H’ contain the components of the continuity 
equation, axial, radial, and circumferential momentum equations, and the energy equation: 



P 

pu 

pv 

pw 

E t 


(2.37) 
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pu 


pu 2 + p + il. 

KC 


F = 


puv + 


'Erz 

Re 


puw + 
r Re 


(E t + P] U — r — 

1 1 ReP^Yo - l) dz 


pv 


puv + 


Re 


G’ = 


pv 2 + ^ 


pvw + -^S- 
K Re 


(E t + P}v - 


k dT 
RePi(Yo-l) 9r 


( 2 . 38 ) 


( 2 . 39 ) 
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H' = 


pv 

r 

P uv P M : B - %z 

r + r Rc 


p 4 - *£+ + - 


Yo 
2pvw 


Re 


t 98 
r Re 


r Re 


| Et + P ) _ l 1 k ?L £ 

» It r RePr ( Yo .i) ^ RcPr(y 0 -l)a 


Rs 


Sr 


_ 5 


_^C : ? +i El 

2 RcPt(Yo- l)r dr dz 


RePr(y 0 .l) 

Lastly, the shear stress terms for an axisymmetric case [21] are: 


Tn = “H 


2— - 1 V-v | 
\ dr 3 


Tee = -pl2^-^ V-v 


k,. 



Ldu 


^ZZ “ 

-P 

2 — 




dz 


where 




Xn = 


-P 

(' i 

ij 

[dw 1 


\dz / 

— 11 

du dv^ 
" + " " 


dr dz 


_ _ id. . du 

V v — (rv) + — 

r dr dz 


For the axisymmetric case, the components of the Current Law, Eq. (2.34), 

dB e 
dz 

dB r dB z 

dz dr 


Jr = “ 
je = 


j, = - f — (rB a ) . 

dr 


(2.40) 


(2.41) 


are: 
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so 


In addition, for the the axisymmetric case B r , B z are identically equal to zero [10], 


jr 

Je 

h 



(2.42) 


Using Eq. (2.42), the components of Ohm’s Law, Eq. (2.30), become 


jr = crEj 

j 0 = Ee = 0 (2.43) 

jz = o E z . 

Using Eq. (2.43), the components of the Induction Law, Eq. (2.33), for the 
axisymmetric case become 

(v x e)t = 0 

(VxE)e = ^ ^ (2.44) 

dz dr 

(V x e) z = 0 

Finally, Eqs. (2.42), (2.43), and (2.44) are combined to obtain the equation 
governing the electromagnetic field: 



Thus, the governing fluid mechanics equations, Eq. (2.36) are uncoupled from 
the governing electromagnetic field equation, Eq. (2.45). The link between the two sets 
of equations is the electrical conductivity, o, which is a strong function of pressure and 
temperature. Once the pressure and temperature distribution are known a can be 
determined using tabular data, empirical models, etc. Given a and appropriate 
boundary conditions for B e , Eq. (2.45) can be solved. The components of j are then 
solved for using Eqs. (2.42)-(2.44). The source terms containing the components of j 
and B 0 appearing in H' can then be calculated. 
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An iterative procedure is thus established between the fluid mechanics and the 
electromagnetic field portions of the flowfield. Both sets of governing equations are 
solved repeatedly until a steady-state solution is achieved. 

2.5 Transformed Governing Equations 

We will now consider the transformation from the physical domain (z,r) to the 
computational domain (^,T|) where 


% = &’*) 
and 

= fl(z,r) 

Applying the chain rule of partial differentiation yields 

d d ^ dr| d 

dz dz dz dq 

a_ _ a_ an d_ 

dr dr dr dt] 

Applying this transformation to the governing equations, Eqs. (2.36), (see 
Appendix B) yields 



d$ dn 


(2.46) 


where the vectors Q, F, G\ H 1 are the Q, F, G\ if vectors in the transformed 
coordinate system. Splitting the F, G\ and H 7 vectors into inviscid (F, G, H) and viscous 
(F v , G v , H v ) components yields 


3Q 0F dF v dG 

— + — + + — + 

dt d£ d£ dq 


^ + H + Hv=0 

dq 


(2.47) 


An assumption has already been made the flow is axially symmetric, causing 
gradients of terms in the circumferential direction to be eliminated. However, the resulting 
axisymmetric Navier-Stokes equations do contain a circumferential or swirl velocity 
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component as well as axial and radial velocity components. These equations can be said to 
model a strongly swirling flow since only variations in the circumferential direction are 
neglected and all other terms are retained in the momentum equations. 

In contrast, a weakly swirling axisymmetric flow can be adequately represented by 
applying the standard boundary layer approximations [12] 


d_ 

dr 


v 0 ,v z *0(1) 

« Op-) where e « 1 


I- 0 * 1 * 


to the axisymmetric Navier-Stokes. However, the radial momentum equation, in the 
physical coordinate system, reduces to 


P v 9 _ dP 
r dr 

a balance between centripetal and pressure forces. In addition, all viscous terms in the 
momentum equations containing derivatives in the strcamwise direction are found to be 
negligible. These are the governing equations used in [1, 10]. 

One of the aims of this study is to determine the effect of swirling flow on the 
current and velocity profiles. Obviously the radial momentum equation is simplified so 
immensely by a boundary layer assumption, no shear stress terms for example, that the 
effect of a large circumferential (or swirl, or tangential) velocity component cannot be 
adequately represented. Instead an assumption that the strcamwise diffusion terms are small 
compared to the normal diffusion terms will be made, i.e., 

d d 

This is referred to as the Navier-Stokes thin layer approximation (TLA) [27]. Using this 
approximation only the viscous terms in the strcamwise direction are eliminated and all 
other terms in the momentum equations are retained. 

The resulting TLA equations contain more information about the flow than the 
boundary layer type equations which are valid only for weakly swirling flows. Thus, the 
TLA will be able to represent a stronger swirling flow than the boundary layer type 
equations yet only neglect the strcamwise derivatives in the shear stress terms. 
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Applying the TLA to Eq. (2.47) yields 


i 9Q 9F 9G 9 G v g fv rv 
± — + — + — + — - + H + H v =0 

J at 9n an 


(2.48) 


where the components of the vectors are 


Q = 


P 

pu 

pv 

pw 

Et 


(2.49) 


F = d A 



o. 


puu' 

+ 

fr' p 

pvu' 

+ 

£- p 


pwu' 



E t + P| u' 


(2.50) 


G « dfi 


pv' 


puv’ + 

ifc r S P 

pvv' + 


pwv' 

(E, + P) 

v’ 


(2.51) 


where u' and v' are the rotated velocity components normal to surfaces of constant £ and n> 
respectively (see Appendix B). 
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Lastly, using the chain rule of partial differentiation the coordinate transformation is 
applied to the equation governing the electromagnetic field, Eq. (2.45) to yield 


+ n, ln(a( , ' rBn + ^ B5 + £)) + k^(o(’ 1,B ’> + 5rB5 + ^| =0 


where 


where 



etc. 


(2.55) 


The relationship between the Jacobian, J, metrics, £ r , £ z , Th, rj z , and inverse 
metrics, z^, r^, z^, r^, of the coordinate transformation as well as details of the preceding 
derivation are contained in Appendix B. 
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Chapter 3 


Numerical Method 


3.1 Introduction 


This chapter discusses MacCormack's [20] implicit, flux-split, Gauss-Seidel line solver 
algorithm. The algorithm is used to solve the finite volume approximation to the governing 
fluid mechanics equations. The implicit and explicit forms of the inviscid, viscous, and 
source term portions of the governing equations are developed. The boundary conditions 
and procedures needed to model solid wall and centerline conditions are discussed and 
presented. An appropriate boundary condition procedure for a subsonic, swirling, viscous 
flow in an axisymmetric nozzle is also developed for the nozzle inlet plane. 

3.2 Implicit Formulation 


The flux-split finite volume approximation to the governing fluid mechanics 
equations in transformed coordinates (recalling Eq. (2.48)) 


+ + 36 + ?e. + ff=0 

^ 3t dq dq 


(3.1) 


will be developed now. 

Since the inviscid flux vectors, F, G, are homogeneous functions of degree one 
[28] we can write 


F=^Q=AQ and G = ^Q=BQ 

dQ 3Q 

where A,B are the rotated Jacobian matrices of F, G and are derived in Appendix C. A 
consequence of this is that the inviscid flux vectors can be split into subvectors. One 
subvector is associated with all positive eigenvalues and the other with all negative 
eigenvalues. The subvectors are then differenced with an appropriate one-sided scheme. 

The flux vectors are split after the Jacobians A, B have been diagonalized in the 
following manner, such that and Aj, are diagonal matrices composed of the eigenvalues 


of A and B, respectively 


A = S'^C^Aa C a R A S 
B = S 1 Rb 1 Cb 1 A b C b Rb S 


where 


A = 


Xi 


X 2 0 
X 3 

0 X 4 


X5 


The matrices diagonalizing A and B, i.e., S, R, and C and the matrix of eigenvalues, A, are 
developed in Appendix C. 

The individual eigenvalues are split by writing 

Xi = Xi + Xi 


where 

1 + Xj + Xi 

and 

1 : 

Xj - Xi 

Xi - 2 

A i 

2 

So that if 

Xi £ 0, 

then X* = Xi 

and 

O 

II 

UT 


Xi < 0, 

then Xf = 0 

and 

>> 

N* 1 

II 

>> 

Using these results we can write 






A = A + + A* 




where the diagonal matrices A + and A~ contain the elements X* and Xi respectively. For 
subsonic flows, using the eigenvalues obtained in Appendix C, we have 


Aa = diag( d A u\ d A u\ d A u\ d A (u’ + c), 0 ) 
Aa = diag( 0, 0, 0, 0, dA(u' - c) ) 
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and for supersonic flows 


A a = diag( d A u', d A u', d A u’, d A (u’ + c), d A (u' - c) ) 
Aa = diag( 0, 0, 0, 0, 0 ) 

Finally, the flux vectors can be written 

F = AQ = S* 1 R a C a A a C a R a S Q 
= S * 1 R a c a (aX + A a ) C a R a SQ 
= (A + + A ) Q 
= F + F 


where 

A + = S* 1 R a C a A! C a R a S A* = S 1 R a ! Q 1 A a C a R a S 

— * /N. -* 

F = A + Q F = AQ 

and 

A = A + + A' 

The G flux vector can be written in the same manner using the eigenvalues and matrices 
diagonalizing the B matrix. 

The finite volume approximation to Eq. (3.1) can now be developed. Implicitly 
differencing the equation using first order accurate differences in space and time yields 


Q n+1 - Q n 


~ C:n+1 _ ^n+1 

+ At I P F + P G 


DGv +1 




Ari 


Ari J 


+ H rn+1 = 0 


(3.2) 


The direction of the difference operators, D, will be discussed as each term is examined. 

3.2.1 Implicit Formulation* Inviscid Flux Vectors 


The finite difference representation of 


+i 


derived by central differencing the expression such that 


(and similarly of 


r»f* n+1 

— ) will be 
At| 


df "* 1 _ FT+iftj - FT-tkj 
A4 A£ “ 


(3.3) 


Thus the flux in the % direction crossing the surface of a control volume located at (ij) 
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enters and leaves through the surfaces located at (i ± 1/2, j). This is illustrated in Fig. 3.1. 
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Figure 3.1 Direction of streamwise flux travel 


The flux crossing the surface located at (i + 1/2, j) in the positive % direction originates at ij 
and the flux traveling in the negative £ direction originates at (i + 1, j). Accordingly, the 
inviscid flux across the surface located at (i + 1/2, j) is 



whereas the flux crossing the surface located at (i - 1/2, j) is 



(3.4) 


(3.5) 


Note that the Jacobian coefficient matrices arc evaluated at the cell edge. Clearly, the- 
net flux crossing the surface of a control volume is proportional to the difference between 
Eqs. Eqs. (3.4) and (3.5), however, the grid point at which each A" is evaluated is not 
explicit because only information at the cell center is available, not at the surface of the 
control volume where the Jacobians need to be evaluated. A method to resolve this 
ambiguity was proposed by Steger and Wanning [28]. For a flux traveling in the positive £ 
direction the Jacobian is evaluated at the point upstream of the surface and for a flux 
traveling in the negative £ direction the Jacobian is evaluated at the point downstream of the 
surface. Thus, the Jacobians would be evaluated at 

A +n v A+ n A-n . A-n 

i+l/2,j A ij A i+1 j 


A -n 


A" n - 


4+n . A +n 

A i-l/2.j A i-lj 
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An alternate method has been proposed by MacCormack [20] and discussed in 
more detail by Candler in [29]. In this method both Jacobians, A n and A +n , are always 
evaluated at the same point The point chosen alternates the half integer indices i+1/2 (i-1/2) 
between i and i+1 (i-1 and i). Equivalently the half integer points at i + 1/2 (i -1/2) can be 
evaluated using the average of the flow variables from the points i and i+1 (i and i-1). 

MacCormack and Candler [30] demonstrated the Stcger-Warming form of flux 
splitting results in excessive numerical diffusion in the boundary layer while the 
MacCormack flux split procedure prevented excessive numerical diffusion. Flows 
containing shocks must be evaluated with a different procedure [31], The MacCormack 
procedure, evaluating the half integer points using the average of the flow variables, was 
used in the present study. 

The inviscid flux vectors are made second order accurate [32,33] by linearly 
extrapolating the value of Q n . Thus the value of value of Q n approaching the i+1/2 surface 

■ * n -*n - *n 

from the upstream (downstream) direction is 3/2Qjj- l/2Q;.i j (3/2Qi + i,j- l/2Qj + 2.j)- 
Pressure switches are used to degrade the accuracy to first order in regions of large 
pressure gradient. For example Eq. (3.4) becomes 


- Aftnj (Q”j + sWiocj (q7j - q7-,.j)) 

+ Aj iaj -f SW i ,3 ftj (Q’’„. J -QU J )j 


with 


SW i+3/ 2 ,j = Max 0.0, 0.5 - 


|Pj+2 - Pj+ll 


Min(Pi + 2 , Pj+i) 


By evaluating A at time level n the inviscid flux can be linearized to yield 

F" +1 = F" + A n (Q n+1 - Q n ) + 0(At 2 ) (3.6) 


In summary, the £ directed inviscid flux at time level n+1 crossing the surface 
(i+1/2, j) is approximated using Eqs. (3.4) and (3.6) as 


Sj 1 + A s:.'j 

= Aj"i/2j Qij + Aj+i/2,j Qi.ij + Aj^i/^j 6Qij + Aj" .,/2j SQ’i.ij 
where SQ = Qi,j - Qij = 


(3.7) 


Similarly the flux crossing the (i-1/2 j) surface is expressed as 
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(3.8) 


Fi-V/2,j= A'" 1/2 j Qi,j + Ajti/2.j QmJ + A £l/2,j &Q?.j + A M/2.j SQi-l.j 

Finally substituting Eqs. (3.7) and (3.8) into (3.3) and gathering terms yields 


^ - -L j[D- A + D. Af ln J (fiQij + ©.j)) (3.9) 

A S A S 

where D„ and D+ are backwards and forward difference operators, respectively. Similarly 
in the rj direction we have 


DG n+1 

Ari 



+ D, B* 1( J («©, + ©j)} 


3.2.2 Implicit Formulation- Viscous Flux Vector 


(3.10) 


D G 

The finite difference representation of — — - — will be expressed at time level n by 

All 

defining 

er 1 = e: + sg! 

can be partitioned into: 1) the terms common to both the cartesian and cylindrical 
formulations and 2) the additional terms arising from the cylindrical axisymmetric 
formulation 


Gv Gvcom Gvcyl 
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J |P 4 + p2V n ) 
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J p3 w n 
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By employing the vector of primitive variables, V 


V = (p u v w T) T 


we can write 


where 


8G V — SGvcom "*■ SGvcyl 

= M.^-8V + M b 8Q 
dx\ 


8V = (8p 8u 8v 8w 8 t| T 


Matrices M, and M b are given in Appendix D. We can change variables from 8V to 

- d V — - 

8Q by defining the Jacobian N where N = — — such that 8V = N 8Q. This yields 

dQ 


8G V = M, - + M b 8Q 

dr\ 


Matrix N is also given in Appendix D. 

Thus the viscous flux in the T| direction becomes 



The middle and final terms in Eq. (3.1 1) are central differenced [34] with the final term 
expressed as 

2- (M b Sfyj = [(M b 8 q)l 0+1/2 “ (Mb SQjio-m] / Ati 

At] 

= [M b i,j+l/2 (5Qi,j+i + SQi.j^ + Mbi,j.i /2 (8Qij + 8Qy.i|/2] / Arj 

In other works, [29,30], the Thin Layer Approximation (TLA, see Section 2.5) 
was also applied during the formulation of the implicit algorithm. Consequently, these 
formulations also neglect the viscous flux in the £ direction in the implicit portion of the 
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formulation. However, the viscous derivatives in the ^ direction could be retained in the 
explicit portion of the foimulation, thus influencing the converged solution [20]. 

✓V 

3F V 

If the £ direction viscous terms, — — , were retained they would be approximated as 

K 

the Gv term was 


& = BIT!. . jl/k + 5?;) 

H A5 A5 l ' 

The TLA is now applied to the both the explicit and implicit derivatives, the implication 
being that gradients of the viscous terms are much larger in the rj direction than in the t, 
direction 

£-(Gv + 5Gv) » iMF? + 8F?) 

AT]' ’ AS, ' ’ 

resulting in 

J^-(Fv + 8F?) = 0 (3.12) 

AO 1 

3.2.3 Implicit Formulation* Source Term 


The source term can be split into inviscid and"viscous” portions. The inviscid 
portion of the source term is due to the additional convection terms that arise due to the 
axisymmetric geometry. The "viscous" portion of the source term contains terms arising 
from the shear stress terms in an axisymmetric geometry as well as terms due to the current 
flow and radiation. Consequently, the source term can be linearized, as the inviscid flux 
terms were, to be 





H"j + Hvi.j + 8H? tj + 8H 


v «.j 


The viscous source terms will be neglected in the implicit side of the formulation, 

8H V i,j = 0, but will be retained in the explicit portion of the formulation, i.e., H V i.j * 0. 
Hence, the physics of the flow will be retained in the explicit portion of the solution (all 
terms are retained in the explicit source term) and therefore influence the converged 
solution. Finally, using these assumptions the source term can be written as 

H”‘ = H"j + CiijSQij (3.13) 
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where 


H'"j = (H + Hv)T.j and 

The matrix can be found in Appendix D. 

3.2.4 Implicit Formulation- Finite Volume Approximation 

Substituting Eqs. (3.9)-(3.13) into Eq. (3.2), yields the finite volume difference 
equation approximating the governing equations 


r _ 3H 
Wi,j ” 

BQ 


i + -Al 


D- A+n . Df A-n , D- p+n , D+ D-n 

A i+l/2,j + 77 A i-l/2 f j + + 


v ij A£ Atj Atj 

1 


+ K,.j + cji.il, 

Anl An / An J / 


5Q" 0 - AQy (3.14) 


where 


Vy = cell volume = 1/J 


AQy - - ^ 
v u 


D- A+n . Df A-n . D- i>+n 

. t A i+l/2,j + A i-l/2.j + "ijn 
Aq Aq Ar| 


D- T>+n . Dj n-n 
“ ^ij+1/2 + — 

Ar| Atj 


Qlj 


Al 

Vij 


-&-G n vi.j + Hy + Hvi.j 
Ar| 


In the present work the grid point spacing in the computational plane is uniform 
yielding A£ = Atj =1. Finally, Eq. (3.14) can be written in a form used in the Gauss- 
Seidel procedure as 


By SQy., + Ay SQ°, + Cy 8®.,., - RHSy 


(3.15) 


where 


By = Al|b3. w * N"y„ + 
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Ay = 

_ AL 

Vy 



+ 


Al_ cs. . 
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▼ II 





+ M 


n 

ai,j-l/2 



Mby-mj 


RHSjj = - D,j 5S + i.j - Eij 5 Q",j + AQ", 


Bij « Km.j 

Vij 


3*3 Boundary Conditions 


Eio 


= __&L At" «*■* . 

v .. a i-1/2,j 

▼ it 


The block-tridiagonal set of equations given in Eq. (3.15) can be written in matrix 

form as 


B JMAX-1 A JMAX-1 CjMAX-1 

• » • 

• • * 



SQjmax 

SQjmax-i 


RHSjmax 

RHS jmaX-1 

Bj Aj Cj 



5 Qj 

= 

RHSj 

b 2 

a 2 c 2 

i 

8Q2 

6Q1 _ 

i 

rhs 3 

rhs 2 


(3.16) 


Boundary conditions representing a solid wall or centerline symmetry conditions on 
the upper and lower boundaries are imposed through the j= 2 and JMAX-1 rows of the 
block matrices shown above. The appropriate inviscid and viscous boundary conditions are 
implemented differently in the explicit and implicit portions of the block matrices and will 
be discussed separately. 
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The boundary conditions imposed at the inlet and outlet planes of the thruster, i= 2 
and i= IMAX-1 respectively, are implemented using a Method of Characteristics (MOC) 
procedure and are also discussed subsequently. 

3.3.1 Boundary Conditions- Upper and Lower Boundaries 

Several different methods could be used to implement the boundary conditions on 
the upper and lower boundaries. One method is to use a fictitious element outside of the 
computational domain. The wall volume element common to the interior elements and the 
fictitious elements represents the physical boundary of the object being modeled, e.g. 
nozzle wall, centerline, etc. Using the boundary conditions and flow information from the 
interior elements the flow conditions in the fictitious elements can be set to yield the proper 
boundary conditions, e.g., no-slip along a solid wall. The flux crossing the surface is 
calculated in the same manner as for all other cells. The resulting flux should represent the 
desired boundary conditions. 

Another method, implemented in the present work, calculates the flux by directly 
applying the boundary conditions to the affected fluxes, e.g., no normal flow across a 
surface representing a solid wall or centerline. Thus, fictitious elements are not used. 
However, to allow for the possibility of investigating several methods of applying 
boundary conditions the notation used in this study also allows fictitious boundary 
elements. The subscripts i j represent a cell located at the ith element in the streamwise and 
-the jth element in the cross-stream direction. Therefore, the fictitious elements not used in 
the present study are elements located at: (i= lwIMAX, j* 1), 

(i= 1=^IMAX, j= JMAX), (i= 1, j= 1=>JMAX), (i=IMAX,j= 1=>JMAX). Thus, in 
a grid described as 60 x 40, IMAX= 60, JMAX= 40, and the computational domain 
actually contains 58 x 38 elements. 

3.3. 1.1 Upper and Lower Boundary- Explicit Boundary Conditions 

The explicit inviscid boundary conditions across a surface of 11= cst will be 
examined first. There is no flux of mass, momentum, or energy across a surface 
representing an impermeable wall or centerline and the normal velocity to the surface is 
zero. These boundary conditions are expressed in the explicit inviscid flux vector as (which 

is obtained by setting v'= 0 in G) 

G = ^ 0 -r^P z^P 0 0 j 1 (3.17) 

The explicit viscous flux boundary conditions are set by the following conditions 
along the lower and upper boundaries 
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+ u'i.2 centerline symmetry 

- u'i,2 solid wall, no slip 

+ v\2 centerline symmetry 

- v'j 2 solid wall, no slip 

( 3 . 18 ) 

+ Wj,2 centerline symmetry 

Wi.l = { 

— Wj t 2 solid wall, no slip 

adiabatic wall/centerline symmetry 

isothermal or prescribed 
wall temperature 

Centerline symmetry velocity conditions state the axial, radial and swirl velocity 
components of the two volume elements on either side of the centerline have the same 
magnitude and direction. An adiabatic wall condition is expressed by a temperature 
equality. There will be no heat transfer across the surface if the temperatures of the two 
cells adjacent to the boundary have the same temperature. The isothermal or prescribed wall 
temperature condition is derived using an arithmetic average of the temperatures of the two 
volume elements adjacent to the wall, T i>wa n = (Tj t i + T^) / 2 . 

For example, the explicit inviscid and viscous boundary conditions for the lower 
surface are implemented through the term RHS^. where 

RHSi,2 = - D it2 5 Q i+u - E i>2 5 Qi-i ,2 + AQi, 2 
by modifying the terms contained in the explicit term AQi,2 




AQi.2 = - At Aj+i/2,2 + ^7 A ’i-l/2.2 + ~ B T.2+l/2 + —■ B i, 1+1/2 

A^ A^ All Aq 



-At 


^Gvi.z + Hi. 2 

An 


( 3 . 19 ) 
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Explicit Inviscid Boundary Conditions 

The explicit inviscid boundary conditions given by Eq. (3.17) are implemented by 
modifying the inviscid terms at the lower boundary, in the T] direction, in Eq. (3.19) given 
by 


— Bu,i/2 + ^ b;.,. 1/2 q u 
\Ari Ati ) 

= ( B T,2+l/2 Qi.2 + ^.2+1/2 Qi.3) - ( B T.l+l/2 Qi.l + ^1+1/2 Qi.2) 


The second bracketed term in the last expression represents the flux crossing the surface 
located at j= 1 + 1/2. Using Eq. (3.17) the resulting flux is 



+ ^.1+1/2 Qi.2 = Gi.1+1/2 = 


0 -r^Pi .2 z^P if 2 0 


0 


TT 


and the metrics are evaluated at the surface j= 1 + 1/2. If the surface is a centerline, the 
metric r^ will equal zero. 

Thus the explicit inviscid flux terms in the ri direction at j=2 become 


D. n+ 


BT.2,1/2 + — B u , 


1/2 


Qi.2 


(At) At| 

= ( B r.2+l/2 Qi.2 + Blj+wQu) - |Gi.l+l/2j 


Following an analogous procedure the explicit inviscid flux at the upper boundary 
is set by modifying the inviscid flux terms at j=JMAX-l 

— B UMAX-l/2 + B i,JMAX-3/2 1 QiJMAX-1 
\Ati Ati } 

= ( B UMAX-l/2 Qi.JMAX-1 + B i.JM AX-1/2 Oi.JMAx) 

- ( B UM AX-3/2 Qi.JM AX-2 + S’i.JM AX-3/2 OiJM AX-1 ) 

The first bracketed term in the above expression represents the flux crossing the j= 
J MAX- 1/2 surface. Using Eq. (3.17) the explicit inviscid flux terms in the T} direction at 
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j=JMAX-l become 


. ~ B UMAX-l/2 + — B i,JMAX-3/2 QiJMAX-1 
IAT| Aq / 


= (Gi.JMAX-1/2) - ( B tjM AX-3/2 Qi.JM AX-2 + B 'i.JMAX-3/2 QiJMAX-l) 


where 


Gi,JMAX-l/2 


0 -r^PijMAX-l Pi.JMAX-l 0 


0 


T 


and the metrics are evaluated at the surface j= JMAX - 1/2. 


Explicit Viscous Boundary Conditions 

The explicit viscous boundary conditions given by Eq. (3.18) are implemented by 
modifying the viscous terms at the lower boundary, in the rj direction, in Eq. (3.19) given 
by 


J2-G., 2 

Ati 



Aq 


+ ^(° vcy1 ^ 2 

Aq 


where 


= {M ai .2 + l/2(Vi,3 - V, 2 ) - Mm.i+i/ 2 (V it 2 - V U )) 
+ (Gvcyl i,2+l/2 — Gvcyl i t l+l/2} 



Gvcyl ~ 


rRe 


0 


-!■ 

2 . 

3 


z $ v 


w 

0 


(3.20) 


The explicit viscous boundary conditions are implemented by modifying the two 
bracketed expressions in Eq. (3.20). The first bracketed expression contains the terms 
common to either a cartesian or cylindrical coordinate formulation and will be modified 
first. 
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Consider the boundary conditions for the lower surface located at j= 1 + 1/2. The 
second term in the first bracketed expression in Eq. (3.20) is 

p2 - Pi 
u 2 -ui 
v 2 - Vi 

W2 * Wi 
_ T 2 ' T 1 _ 

Using the explicit viscous boundary conditions given in Eq. (3.18) this expression can be 
written as 


»i.l+l/2(Vi,2 - V U ) =M,i t i + i/2 




" P2 ' 
u 2 

M ai ,i + i/2(Vi ( 2 ~ V U ) = lvCi.i+i/2 v 2 = ^ai.1+1/2 Vi.2 (3.21) 

W2 

. T 2 _ 
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where 


0 0 0 0 


0 


M a i,l+l/2 


J Re 


o tup 4 tufo o 0 

0 tup 2 tuPl 0 0 

0 0 0 t u p 3 0 


0 0 0 0 tr 


1 


Pr(ito-1) Ps J 


Following an analogous procedure the second bracketed expression in Eq. (3.20), 
containing the additional viscous terms arising from the cylindrical coordinate formulation, 
can be modified to implement the viscous boundary conditions along the lower surface. It 
should be noted that the terms for the surface j + 1/2 are evaluated using the average of the 
cell variables at j+1 and j. 

At the lower surface the boundary conditions are implemented by modifying the 
term at the surface 


✓s 

G 


vcyl i,l+l/2 


(3.22) 


After applying the viscous boundary conditions given in Eq. (3.18), and L'Hospital's Rule 
the last expression can be written as 

r o i 


o 


|j. 

Gvcyl i,l+l/2 = 


7 dv 
3 z *dF 
dw 

o 


(3.23)- 


Substituting Eqs. (3.21) and (3.23) into Eq. (3.20) yields the modified expression 
for the viscous terms on the lower boundary 


^-Gvi,2 = (M.i,2+l/2(Vi,3 - V ii2 ) - Nf,i.i+i/2 Vu) 

At] 


+ Gvcyl i,2+l/2 Gvcyj i.1+1/2 


(3.24) 
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Following an analogous procedure the explicit viscous boundary conditions at the 
upper boundary are set by modifying the viscous flux terms at j=JMAX- 1 


-^-Gvi.JMAX-l 

Ar| 



Mai.JMAX-l/2 — VijAMX-1 1 
All 


+ 


-^-(GvcylUMAX-l 

AT] 



«i.JM AX-1/2 


(ViJMAX “ Vi,jMAX-l) ~ M a i,jMAX-3/2(Vi,JMAX-l 
+ (Gvcyl i,JMAX-l/2 ~ Gvcyl i.JMAX- 3 / 2 ) 


~ VijMAX-2) 


The modified viscous flux terms at j=JMAX-l are 

Gvi,JMAX- 1 = (- N^li.JM-l/2 V i.JM-l - Mai,jM-3/2 (Vi.JM-l “ VijM- 2 )) 

* 

- Gvcyl i.JMAX-3/2 

^ai.JM- 1/2 contains the same terms as M* ; 1+1/2 but is now evaluated at iJMAX-1/2. 


The tj derivative terms appearing in the explicit viscous source term, H v , are evaluated 
using central differences just as was done for the G VC yi terms. 


3. 3. 1.2 Upper and Lower Boundary- Implicit Boundary Conditions 

The implicit inviscid boundary conditions across a surface of ti= cst will be 
examined first. The following conditions are specified to maintain impermeability of a solid 
wall or enforce centerline symmetry at the lower boundary 


Spu 

= 8pi, 2 

8pu’i,i 

= Spu’i.2 

+ 5pv' it 2 

centerline symmetry 

- Spv’i .2 

solid wall 

8pw it i 

= 5pWi,2 


= 8E ti(2 


(3.25) 
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where u’ and v' are the rotated velocity components (See Appendix B) tangential and 
normal to a surface where T|=cst . 

The implicit viscous flux boundary conditions are set by the following conditions 
along the lower boundary 


8u' : 


8v'i,i = 


5w, i = 


5T U = 


1 + 8u'i,2 

centerline symmetry 

\- Su’i.2 

solid wall, no slip 

+ 8v'i,2 

centerline symmetry 

- 8v’i,2 

solid wall, no slip 

r 

+ Swi, 2 

centerline symmetry 

- 8wj,2 

solid wall, no slip 

I +6T U 

adiabatic wall/centerline symmetry 

| -ST ii2 

constant temperature wall 


(3.26) 


These conditions are obtained by applying the operator 8 = At ^ to the explicit viscous 
boundary conditions in Eq. (3.18). 

The implicit boundary conditions for the lower and upper boundaries are 
implemented by eliminating the C u 2 and Bumax-i. matrices, respectively. The terms 
contained within them are incorporated into the A 1,2 and AjjMAX-i matrices using the 
implicit boundary conditions to form the modified matrices A^, AjjMAX-i- 

Implicit Inviscid Boundary Conditions 

The lower surface boundary conditions will be illustrated first The inviscid terms 
in the cross-stream direction at the lower surface are contained in 


Ai. 2 => ^(Bt,2+i/2 “ Bu+itf) 


Ci,2 


y 

^M-B 

'*} 


+n 


v^r Di - i+i /2) 


The inviscid terms at the lower surface, j= 1+1/2, are combined along with the 
boundary conditions and placed into A 2 as 
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a'. 2 =» + Bl" wm )| (3.27) 

Vij 

= v^fc 2 * 1 ' 2 " ( 5B '- 1+1/2 )) 

where SBy+i# implements the implicit inviscid boundary condition 

5Qj,2 = [ 0 -r^5Pi,2 z^8Pj t 2 0 0 J 1 

and is given by 


0 


0 


0 


0 


0 


-r^ a P r^ u P r^ v P r^ w P -r^ P 


881,1+1/2 = 


a P -z^ ^ u -z^ v P -z^ w P z^ P 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


a . » 2 * . v? . ± tt* p - (y- 1) 

All properties are evaluated at j = 2 and the metrics are evaluated at the j = 1 + 1/2 
surface. Thus Eq. (3.27) is the implicit analog to Eq. (3.17). 

In an analogous fashion the inviscid terms contained in Ajmax-i are obtained by 

/S 

placing the inviscid term from Bj,jMAX-i into Ajmax-i and applying the boundary 
conditions to obtain 5Bj,jMAX-i/2 such that 


Ajmax-i = ^(6Bi,JMAX-l/2 - (BijMAX-3/2 )) 

All properties in SBijm AX-1/2 are evaluated at j = JMAX-1 and the metrics are evaluated at 
the j = JMAX - 1/2 surface. 

Implicit Viscous Boundary Conditions 

The lower surface boundary conditions will be illustrated first The viscous terms in 
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the cross-stream direction at the lower surface are contained in 


a 2 =*- ^:|m^2 +1/ 2 n?, 2 + Ki.um N?,2 - 




The viscous terms contained in C 2 at the lower surface, j = 1+1/2, are combined 
along with the boundary conditions and placed into A 2 , resulting in 


A*2=>- y^7 

V U 


^ ai , 2 + 1/2 N ?,2 + < 1 + 1/2 N ?.2 ~ 


fMbi,2+l/2' <.1+1/2 


(3.28) 


M‘ ai , 1 + 1/2 is the same matrix used in the explicit calculation. Mb;,i+i /2 is obtained by 
applying the boundary conditions and LUospital's Rule as was done for the explicit 
viscous boundary condition. 

In an analogous fashion the viscous terms contained in Ajmax-i are obtained by 

«<*S. 

placing the viscous terms from BjjMAX-i into Ajmax-i and applying the boundary 
conditions to obtain 


Ajmax-i => - 


Y^( M ai,JMAX-l/2 N?,JMAX-1 + M*i,JM AX-3/2 N i!jMAX-l 


_ AL 

Vij 


M bi.JMAX-l/2 ‘ Mbi,JM AX-3/2 


AX-1/2 is the same matrix used in the explicit calculation. Mbj,jMAX-i/2 is obtained by 
applying the boundary conditions and is found to be 


M bi.JMAX-l/2 = [°] 


Thus, the Gauss-Seidel expression to be evaluated at j = 2, i.e.. 


B i>2 5Qi,3 +A it2 8Qu+ Cj, 2 8Q itl = RHSi. 2 


can be expressed as 
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Bi.2 5Qi3 + A, 2 5Qi,2 = RHS* 2 

where A *2 and RHS* 2 contain the modified implicit and explicit terms, respectively. The 
modified terms implement the lower surface boundary conditions. 

The Gauss-Seidel expression to be evaluated at j = JMAX-1 

Bi.JM-l 5Qj,JM + Ai t JM-l 5Qi,JM-l + Ci.JM-l 5Qi,JM-2 = RHSi,JM-l 


can be expressed as 


A i, JMAX-1 SQi, JMAX-1 + Cj, JMAX-1 5Qi,JMAX2 = RHS 


i, JMAX-1 


where Ajjmax-i and RHS*j Max .i contain the modified implicit and explicit terms, 
respectively. The modified terms implement the upper surface boundary conditions. 

The final form of the block-tridiagonal system, Eq. (3.16), is 


— * ~ 

A JMAX-1 C JMAX-1 


SQjmax-i 


BHSjmax.j 

Bj Aj Cj 


6Qj 

= 

RHSj 

§ 2 A* 2 _ 

i 

8Q2 

i 

RHSj 


where superscript * indicates the term has been modified to implement the inviscid and 
viscous boundary conditions. 


3.3.2 Boundary Conditions* Inlet and Outlet Planes 

The Method of Characteristics (MOC) procedure about to be discussed was 
developed by Rai and Chaussee [35] as a variant of the method proposed by Chakravarthy 
[36]. Discussions and examples of implementing the following MOC procedure can also be 
found in [37, 38, 39]. 

The five eigenvalues of the E flux in this study are given by (see Appendix C): 

^1,2,3= d A u ’ 

X 4 = d A (u’ + c) 

^5= d A («' - c) 
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These eigenvalues imply waves travel in specific directions (see Table 3.1). For example, 
in a subsonic inlet flow there are 4 positive eigenvalues given by Xj , 2 , 3 , 4 - The fifth 
eigenvalue is negative and conveys information from the interior of the flow, upstream to 
the inlet plane. The equations containing the positive eigenvalues, which would propagate 
information from outside the boundary to the inlet plane (in the positive £ direction), need 
to be discarded and replaced with boundary conditions since the characteristics coming 
from outside the boundary are not defined [36]. 

The MOC procedure will be briefly described. First the characteristics propagating 
information from outside the computational domain to the boundary are eliminated. Then 
the finite difference/volume representation of the governing equations is multiplied by the 
eigenmatrix normal to the boundary. This yields a decoupled set of characteristic equations. 
Multiplying the characteristic equations by a selection matrix allows the characteristic(s) 
within the computational domain to be retained. The eliminated characteristic equations are 
then replaced by appropriate boundary conditions. 


3.3.2, 1 Inlet and Outlet Planes- One- dimensional MOC Example 


The MOC procedure will be illustrated using a simple example [36, 39]. First the 
characteristic equations for the one-dimensional Euler equations will be obtained. The quasi 
one-dimensional Euler equations can be written in vector form as 


where 


dQ dE 
dt + dx 


H 


0 = 

’ pa ' 

pua 

E= 

pua 

(pu 2 + P)a 


ea 


(e + P)ua 


0 


H= 


pda. 

dx 

0 


In these expressions a= cross-sectional area.Defining the Jacobian matrix as A 

, dQ , A dQ _ u 

produces - 5 - + A = H 


dE 

dQ 


Transforming A to a diagonal matrix using the similarity transforms 
A = L ' 1 A L yields 

§ + L ' 1 = H 

dt dx 

where A is the diagonal matrix containing the eigenvalues X 1 , 2 , 3 = u, u+c, u-c and L is a 
matrix composed of rows of eigenvectors. 
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Table 3.1 


Required Boundary Conditions at Inlet 
and Outlet Plane Boundaries 


Speed 


c > u' > 0 


c > u’ > 0 


u' > c 


Specified 

Signal Description Boundary 

Propagation Conditions 



One 

characteristic 

equation 

— \ 


Subsonic Inflow 
Xl. 2 , 3.4 > 0 

Xs < 0 


P t , T t> v/u, w/u 


Four 

characteristic 

equations 



Outflow 

Plane 




Subsonic Outflow 
Xl, 2 , 3,4 > 0 

Xs < 0 


p 


Five T t 
characteristic I 

Outflow 

Plane 


None, 

Supersonic Outflow Extrapolate to 
Xi.2.3,4.5 > 0 obtain boundary 

values 
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Multiplying the above by L and defining 


8Q _ t §Q 

at at 


and H = L H yields 


ag 

at 


a^ = ® 

ax 


which is equivalent to the three decoupled characteristic equations 


( 3 . 30 ) 


aqi 

at 



Each of the three equations has a specific direction and governs wave propagation in that 
direction. Thus there are left and right running waves in the flow. 

Next a flux split algorithm to solve the one-dimensional Euler equations will be 

considered. As developed earlier for MacCormack's algorithm the flux vector, E, in the 
one-dimensional Euler equations can be split into two pans with positive and negative 
eigenvalues 

E = E + + E 


Using the Steger and Warming form of splitting 


and 


A = A + + A 



In addition 


A = A + + A- 


with A + = L' 1 A + L A" — L' 1 A L 

and since E is a homogeneous function of degree one we can write 

E = E + + E = ( A + + A) Q 

By employing Euler implicit differencing in time, the flux split system can be 
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written as 


where 


l?A + 

I - At D + At ^ 5 — 
dx 



-At RHS 


(3.31) 



RHS 


dE + 9E 
dx + dx 



The MOC will now be illustrated for subsonic flow at inlet and outlet boundaries. 
The first step is to neglect the incoming characteristics at the inlet in both the explicit and 
implicit portions of the algorithm. At the inlet the discretized equations reduce to 


I - At D + At 


av 

dx 


5Q = -At j^- - H = -At RHS 


and at the outlet boundary 


I - AtD + At 


dA + 

dx 


5Q = -At 


' 9E* 

dx 



-At RHS 


In the next step the modified governing equations at the inlet are multiplied by the 
eigenmatrix normal to the boundary, L. Using the definitions 


A* = L 1 A L 


and 


SQ = L 6 Q 


yields 


I - AtD + At 


9 A' 
dx 


8 Q - -a.l(|: . h)" 


which is equivalent to the characteristic equation given in Eq. (3.30). Writing the last 
equation in terms of 8 Q yields 


I - AtD + At 


d A 


dx 


_ \ 
8 Q = -AtL HE _ . H 
dx 


(3.32) 


The next step is to multiply the characteristic equations given by Eq. (3.32) by a 
selection matrix, N', such that only the characteristic corresponding to the X 3 eigenvalue is 
retained at the inlet Thus N' is given by 


N- = 


000 
000 
Loo 1 J 
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Similarly, at the outlet in a subsonic flow the selection matrix, N + , must select the two 
characteristics corresponding to the two eigenvalues and 


N + = 


100 
010 
L0 0 0 


Thus the expression at the inlet boundary, Eq. (3.32), becomes 


N'L 


I - At D + At 


3 A 


3x 


8Q = -AtNL - H 

OX I 


(3.33) 


whereas at the outlet boundary 


N + L 


I - At D + At 


3 A + 


3x 


8Q = -AiN + l|^- - hJ 


The final step is to replace the characteristics that were removed with the appropriate 
boundary conditions. For a subsonic inlet flow the stagnation pressure, P t , and stagnation 
temperature, T t , are often specified. 

The boundary conditions are expressed in vector form [36] by 

H(q) h (P t T t 0) T (3.34) 

Performing a Taylor Series expansion of £i with respect to time yields 


But 


“*n+l 

Cl 


. 3Q ~Z n A 3Q3Q 

Cl + At -r — + '•• s fl + At — — — + • 


3t 


3Q _ Q n+1 - Q n _ 8Q n 


3Q 


3t 


3t 


At At 


so 


And hence 


-*n+i — n ao 

Cl = Cl + -^8Q 

3Q 


“*n+l 

8Ci= G - Cl = 8Q 
3Q 
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•n+1 


To accelerate convergence to the desired steady-state solution, Q is held constant 
at the specified inlet distribution of inlet 


— 3Q - — — n 

5Q = ^5Q= Qimc- Q 
3Q 

Finally, Eq. (3.35) is added to Eq. (3.33) to yield 


(3.35) 


^ + N*l(i - At D + At 3^-) 
dQ \ ax | 


_ —n /0F ' \ n 

5Q =£2iniet- O -AtN L - H (3.36) 

ox , 


This completes an example of how the MOC procedure is applied to the one- 
dimensional Euler equations. 


3. 3. 2, 2 MOC Applied to Present Numerical Procedure 


The MOC procedure as applied to the present numerical algorithm will be illustrated 
now for a subsonic inlet flow (see Fig. 3.2). Consider the control volume at i= 2 whose 
left edge is the inlet plane of the computational domain 


Inlet 

Plane 



■ 

i= 1 


i= 2 


i= 3 


j 


Figure 3.2 Inlet plane boundary condition schematic 

The first step of the MOC procedure is to modify the numerical algorithm by 
removing the positive traveling characteristics at the inlet plane. Thus the implicit flux term, 
E 2 .j , and the corresponding explicit flux term contained in RHS 2 ,j are set equal to zero in 
Eq. (3.15) 

Eli = =>° 


56 



In the next step the modified numerical algorithm is multiplied by the eigenmatrix 
normal to the boundary, L, and the selection matrix for a subsonic flow at the inlet plane, 
N*. From the definition of E 

E = AQ = S' 1 R A C A A a C a R a S Q = L^AaLQ 
we have L 1 = S 1 R* 1 C A ! L = C A R A S 

and N- = diag(0 0 0 0 1) 

yielding 

N* L [LHS Modified Eq. (3.15)] = N‘L [RHS Modified Eq. (3.15)] (3.37) 

The final step is to replace the characteristics that were removed with the appropriate 
boundary conditions (since premultiplying the modified governing equations by N~ L 
yields one equation for the five variables Q). For a subsonic swirling inlet flow Dutton [40, 
41] found that specifying the stagnation pressure, P t , stagnation temperature, T t , ratio of 
swirling velocity to axial velocity, <)>= w/u, and the ratio of radial velocity to axial velocity, 
0= v/u, yielded the best results. Chang [39, 42] also specified these boundary conditions 
for his modeling work of swirling flow. 

The boundary conditions are given by 

5[q) = (P, T t a * 0) T (3.38) 


The Jacobian of Eq. (3.38), supplies the remaining four conditions at the inlet. Thus 

BQ 

the term A^"j /2 .j SQl.j is removed and replaced with boundary conditions given by 


^8Q 


dQ. 


. Adding 5Q = ftiniet - ft to Eq. (3.37) results in five equations for the 

idQ jhi dQ 

five variables contained in Q. 

An analogous MOC procedure is applied at the exit plane of the thruster if the flow 
is subsonic. However, the selection and boundary condition matrices are different because 
the number of characteristics traveling to the outlet plane from inside the computational 
plane is different than at the inlet plane. For a subsonic flow there are four characteristics 
conveying information to the outlet plane from the interior of the computational domain. 
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Thus, the selection matrix is given by 


N- = diag(l 1110) 

The one characteristic moving in a negative direction is replaced by specifying the 
static pressure, P. The single boundary condition is given by 


= (0 0 0 0 P) T = 0 


If the flow is supersonic at the outlet plane all of the five characteristics are 
conveying information to the outlet plane from the interior of the computational domain. 
Thus values of 8Q at the outlet plane are obtained by extrapolation from interior values. 
However, in the present study extrapolation is used at the exit plane in both subsonic and 
supersonic regions. 

3.4 Solution Procedure 

To implement the Gauss-Seidel line relaxation procedure the explicit solution, 
AQ" j, is first obtained at all points in the flowfield using data from the current time level, n. 
The flowfield is then swept in the axial direction from the nozzle exit plane to the inlet and 
then swept in the opposite direction. At each axial location, i, the block-tridiagonal matrix 
equation given by Eq. (3.29) is solved at all j using a routine described in [27], The most 

current information available for 8 Qi.T 1 or SQi+i 1 is used during the sweeps. After 
sweeping the field in both directions the solution is updated using 

Ql*‘ = Q" j + 6Q" (3.39) 


The procedure is continued until the solution has converged, i.e., 8Q n * 0. 


It has long been an axiom of mine that the little things 
are infinitely the most important. 

Sir Arthur Conan Doyle, 1925 
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Chapter 4 


Grid Generation 


4.1 Introduction 

Complicated two- or three-dimensional boundary shapes present several challenges 
when attempting to obtain a finite difference solution for a fluid mechanics problem. For 
instance, irregular boundaries will not coincide with a regularly spaced mesh as illustrated 
in Fig. 4.1. 



Figure 4.1 Mismatch between an irregular boundary and a 
rectangular grid 


Hence, it would be necessary to interpolate boundary conditions for grid points not falling 
on a grid point Because boundary conditions govern the solution obtained in the interior of 
the region, they must be accurately applied and interpolation would produce errors which 
will be propagated throughout the interior region. 

It may be desirable to "pack" grid in a particular region of interest, e.g. near a wall 
to resolve the boundary layer flow as illustrated in Fig. 4.2. 



Figure 4.2 Grid clustering near a solid surface 
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These difficulties can be overcome by using the coordinate transformation, 
discussed in Section 2.5, to transform the physical plane into a computational plane as 
illustrated in Fig. 4.3. 


Y 


T\ 



$ 


Physical Plane 


Computational Plane 


Figure 4.3 Coordinate transformation between physical plane 
and computational plane 


The coordinate transformation results in an irregularly shaped object in the physical plane 
becoming a rectangular object in the computational plane. Because the flowfield boundary 
is now rectangular, boundary conditions may be applied without interpolation between grid 
points which in turn allows the clustering of grid points in any region in the physical plane. 

Another possible source of error is the skewness of a cell in the physical plane, 
especially near the boundaries. By providing the ability to compare grids with a variety of 
cell skewness, especially grids with varying degrees of boundary grid orthogonality, the 
effect of the grid on the numerical solution can be determined. 

The ability to generate a numerical grid for a wide variety of converging-diverging 
nozzles was needed to accomplish the objectives of this work. Such a scheme should 
provide for orthogonal grid near boundaries and the ability to "pack" grid in regions of 
large flow gradients. A grid generation scheme possessing these properties is developed in 
this chapter by solving a set of elliptic partial differential equations. The initial scheme is 
found to maintain the boundary distribution of grid points throughout the interior region if 
Poisson's equation is used in place of Laplace's equation. However, the resulting grid is 
not orthogonal to the boundaries. The source terms, or control functions, appearing in the 
Poisson equations are then modified to provide boundary orthogonality and spacing control 
of grid points near the boundary. 

4.2 Previous Research 


One method of numerically generating the grid coordinates involves the solution of 
the following Poisson equations 
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^xx + 4yy — P(£**l) 
TUX + T|yy = Q'(^Tl) 


By switching the independent and dependent variables the equations are transformed to £,r| 
coordinates resulting in 


a - 2 p + y x nT1 = -J 2 (P’(SJ1) + Q’(^ti) x^) 

ay^ - 2py^ + Yynn = -J 2 (p'(5/n) y$ + Q’(t'n) y^) 

a = x 2 + 

P = x^x n + y^T! 

T = x| + y| 


(4.2) 


where J denotes the Jacobian of the coordinate transformation 


J = ffqj = ^ ^ (4 ' 3) 

Thompson.et al. [43] employed this method using various exponentional functions 
containing adjustable parameters for the source terms P’ and Q\ These parameters allowed 
the interior grid spacing to be controlled to some degree. However, as pointed out by 
Thomas and Middlecoff [44], the proper choice of parameters was geometry dependent and 
required experimentation. 

The spacing of grid points along the physical boundary can be set as desired. While 
this spacing provides Dirichlet boundary conditions for the physical coordinates in the 
transformed plane (Eqs. 4.2) it is difficult to maintain this spacing in the interior of the 
region. When P’ and Q' = 0, i.e., when Laplace equations are used to solve for the interior 
grid distribution, internal spacing control is lost and the boundary grid point distribution is 
not maintained in the interior. Instead, Poisson equations must be used to maintain the 
boundary grid distribution in the interior by developing expressions for the source terms P’ 
and Q' on the boundary. 

Thomas and Middlecoff [44] choose source terms such that Eqs. (4.2) possessed 
exponential solutions but were not exponential functions themselves and hence had no 
parameters to adjust 


p' = ofenlfe + 

q’ = vte.TiHnx + rg)- ^ Lnli 


(4.4) 
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Eqs. (4.4) are then substituted into Eqs. (4.2) to obtain 


a ( + ♦ - 2 p x$, + y( Xqq + V Xti) = 0 5) 

a ( + <J> y^) - 2 p y^q + y( + V yq) = 0 

Constraints that the transverse coordinate curves be locally straight and orthogonal to the 
boundary were then developed. Along £= £(x,y) = constant (See Fig. 4.4) we can write 

d£ = 4* dx + dy = 0 

or (4.6) 


dy] _ _5* 

J.. _ 


q= const = tj b 

i i i i I 

/ £= const 


Figure 4.4 Grid boundaries in the computational plane 

But from the coordinate transformation (see for instance Anderson et al. [27]) the 
metrics are = y-q/J and £ y = -Xq/J. Thus, the slope can be expressed as 


(&) = 
'dx '^=cst 


yn 


(4.7) 


Similarly, along T|= T](x,y) = constant 

/dy] = Jlx = n 

Idx L=csl Tiv 


(4.8) 


For % and rj surfaces to be perpendicular the product of their slopes must equal -1 
and local orthogonality at a boundary ti= constant can be obtained using Eqs. (4.7) and 
(4.8), i.e., 

I! MS)-- 

or (4.9) 

Xq + y^yri = 0 = p 
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The constraint that in the neighborhood of the boundary q= T| b , as illustrated in 

Fig. 4.4, the transverse coordinate lines constant be locally straight (i.e., have zero 
curvature) by requiring the change in slope be zero 


a_ 



= 0 


(4.10) 


Expressions can now be derived for 0 and y. Eqs. (4.5) can be combined to 
eliminate y between the two equations. Along a boundary tj= we can then apply the 
condition of local orthogonality, Eq. (4.9), and zero curvature, Eq. (4.10). A limiting form 
of Eqs.(4.5) can then be solved directly for the parameter <|> that is valid at the boundary ri= 
t] b , i.e., 

♦ = • h x u + /( x I + 4) (4 - H) 

An analogous procedure can be followed to obtain an expression for y along S=^ B , 
or by merely substituting <J>,^ for y,T|, respectively, in Eq. (4.1 1). Hence, 

V * - ( x ti *titi + y-n y^Ti) /(x* + y^J (4.12) 


The parameters <() and y along the boundaries can now be centrally differenced and 
evaluated from the given grid point distribution established on tj=t] b (horizontal 
boundaries) and £=£ B (vertical boundaries) t respectively. 0 and y in the interior are 
determined by interpolation. The former, 0, by linearly interpolating along vertical lines of 
constant and the latter, y, by linearly interpolating along horizontal lines of T|= constant. 
Once interior values of 0 and y are obtained, Eqs. (4.5) can be solved iteratively for the 
interior values of x,y by SLOR. 

We will now apply this procedure to the physical domain shown in Fig. 4.5 and 
used by Thomas and Middlecoff [44] 



Figure 4.5 Thomas and Middlecoff boundary in the physical 
plane 
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Figure 4.6 Grid distribution for Thomas & Middlecoff 
geometry; equal spacing on boundary BC, 

(a) Laplace solution (b) Poisson solution 
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In this case, the arc BC (the nozzle wall) was generated by the superellipse 
= 1, where n =3, b/a = 1.0 / 1.4. Along segments OC, and AB boundary 

values (x,y) were clustered near the points C and B, respectively. Equally spaced boundary 
values were specified along boundary segments OA, and BC. 

Figure 4.6 shows the results for P’=Q'= 0 (Laplace’s equation), and for P' and 
QV 0 (Poisson's equation, using the results developed above). In the Laplace solution it 
can be seen that the boundary point distribution is not maintained in the interior, while the 
Poisson solution maintains the boundary spacing throughout the interior. 

Good grid orthogonality at the boundaries is obtained on all boundaries except 
segment OA as point A is approached. Non -orthogonality occurs because the boundary 
conditions derived for P' and Q' assumed coordinate lines were locally straight and 
orthogonal at the boundary, but these conditions were not enforced in the numerical 
routine. 

Figure 4.7 shows the results for the Laplace and Poisson solutions when the 
boundary point distribution on segment BC is not equally spaced (all other boundary 
segment spacings remained the same). The boundary point distribution is again maintained 
in the interior in the Poisson solution when compared to the Laplace solution. However, 
orthogonality along boundary segment BC is very poor. Thus, choice of the boundary 
point distribution greatly affects the orthogonality of the solution using the P' and Q’ 
relationships developed above. 

The Poisson procedure that has been described was applied to a geometry 
representative of current low-power arejet thrusters [6]. Figure 4.8 shows the results for 
two different boundary point distributions along the top and bottom boundaries. The first 
result is for the case in which the top and bottom distributions are the same and are 
clustered about the center of the nozzle throat The second is for the same bottom boundary 
distribution but with tighter clustering around the nozzle throat along the top boundary. In 
both cases the left and right boundary distributions are the same and are proportioned to 
resolve the wall boundary layers in viscous flows. 

In both cases the boundary point distribution has been maintained in the interior. 
Orthogonality along parallel surfaces has been maintained but along the sloping walls there 
is poor orthogonality. In fact the "vertical" coordinate lines are essentially vertical 
regardless of the slope of the wall they intersect. Magnified views near the wall and at the 
throat, see Figure 4.9, verify that the grid is not orthogonal at the boundary. It should also 
be pointed out that the solution obtained for different amounts of clustering about the nozzle 
throat does exhibit slightly better orthogonality at the wall. This is by virtue of the varied 
amounts of clustering used on the top and bottom boundaries. 

If a numerical routine is developed that provides orthogonality at a boundary we 
should expect to see vertical grid lines, for a case where the top and bottom boundaries 
have the same point distribution, that have an "S" shape appearance (see Fig. 4.10). 
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Figure 4.7 Grid distribution for Thomas and Middlecoff [44] 
geometry; clustered spacing on boundary BC, 

(a) Laplace solution, (b) Poisson solution 
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(a) 



Figure 4.8 Poisson grid distributions for an arcjet thruster 
geometry with, (a) Same point distributions, (b) 
different point distributions, on top and bottom 
boundaries 
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Figure 4.9 Magnified view of Poisson grid distributions for 
an arcjet thruster geometry ( (1) upstream nozzle 
wall, (2) nozzle throat) with, (a) same point 
distributions, (b) different point distributions, on 
top and bottom boundaries 
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Figure 4.10 Expected orthogonal boundary grid distribution 
for an unclustered boundary grid pt. distribution 

On the other hand if we specify the top boundary point distribution to be clustered 
tighter than the bottom boundary distribution we should expect to see vertical grid lines 
such as those drawn in Fig. 4.1 1. 



Figure 4.11 Expected orthogonal boundary grid distribution 
for a clustered boundary grid pt. distribution 

Thus, to improve orthogonality throughout the interior of the region it is 
advantageous to cluster the top and bottom boundaries differently and eliminate the "S" 
shaped grid that can be expected to be produced by an algorithm ensuring orthogonality at 
the boundary. 

4.2 Wall Orthogonality Modifications 

Because the numerical algorithm is not designed to ensure that the grid is locally 
orthogonal at the boundary another method was needed to improve the approach used. 
Note that P' and Q', as calculated above on the boundaries, do not require any information 
from the interior of the region. Examination of Eqs. (4.1 1) and (4.12) reveals 4> and y only 
require information along their respective boundaries, e.g., <|> is calculated on a line of 
constant T| and is centrally differenced using the surrounding values of x and y in the £ 
direction and y is calculated on a line of constant £ and is centrally differenced using the 
surrounding values of x and y in the Ti direction. The control functions contain no 
information from the interior region and are unable to produce orthogonal grid near a 
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boundary unless a fortuitous choice for the boundary grid point distribution is made. 

Thompson [45, 46] has derived source terms from the same governing equations 
used above but which require interior grid point information. As will be shown these 
source terms produce a grid which is locally orthogonal and maintains the boundary point 
distribution in the interior region. 

Table 4.1 contains a summary of the notation differences between the work by 
Thomas & Middlecoff [44], the current research, and Thompson [45, 46].Using the 
Thompson notation we can rewrite Eq. (4.5) as 

M 2 ( x « + Px $) * 2 f i- f n x ^ + N 2 ( x ™ + Q x ^) = 0 (413) 

M 2 (y^ + p y$) - 2f l • riiy^ + l r d 2 (y™ + Q *i) = 0 

Similarly the orthogonality boundary condition, Eq. (4.9) can be written as 

r^ • rtj = 0 (4.14) 

Multiplying the first of Eqs. (4. 1 3) by i and the second by j and adding yields 

| r n| 2 (%5 + p %) - 2r^ • r T1 r^ T1 + |r^| 2 (r^ + Qr^) = 0 (4.15) 

The orthogonality condition, Eq. (4.14), is imposed next on Eq. (4.15) to obtain an 
equation for P and Q valid on the boundary 

M 2 (%S + pf d + N 2 ( f ’Pi + Q f n) = ° < 41 « 

■* 

We now find — — and — — of Eq. (4.16) and obtain the following expressions 

hi 2 N 2 

for P and Q 

P = . • Si . h : hi = 

W 2 M 2 

o _ ^ M = 

N 2 ' W* 


♦ . 


r 5 * r nn 


(4.17) 


V - 




r d 2 


(4.18) 
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Table 4.1 


Grid Generation Notation (Rosetta Stone) 


Thomas & Middlecoff 

Current Research 

Thompson 


P’ 


J 2 


J 2 


Q’ 

<3(e.Tl)|r 6 | 2 

J 2 


J 2 

a 

x ^ + 

M 2 

Y 

x i + A 

w * 

P 

x^ + y&n 

h ■ ? -n 


J 

x^ti - x^ 

- (X 5 x„ + y 5 y w ) 

<t» 

- • hs 

N + 4) 


i*i 2 

* ( x ti + Yn yTn) 
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•» 
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+ $ 
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Thus the control functions derived by Thompson, P and Q, contain the control 
functions used in the Thomas and Middlecoff approach, <{> and \j/, as well as additional 
terms. These additional terms are related to the radius of curvature of the boundary [47] and 
are evaluated using the interior point distribution. 

For example, in P the second term is 

h • = X S + yy^n 

M 2 ( x n + $ 

P is evaluated along two boundary lines of constant r\ corresponding to the outer wall, and 
the cathode surface/centerline of the arcjet. x varies along this line of constant T|. Hence, as 
discussed previously for <|> and \j/, the derivatives x^ and y^ can be directly evaluated from 
the given boundary point distribution using central differences. The derivatives in the r\ 
direction are evaluated along a line of constant £ using one-sided differences into the 
interior of the region. 

The grid generation algorithm was modified to include the source terms shown in 
Eqs. (4.17) and (4.18). P and Q are calculated on the boundary using an initial interior 
point distribution based on the boundary point distribution. Interpolation from the 
boundary values of P and Q are used to obtain P and Q in the interior. Equations (4.13) (as 
were Eqs. (4.5)) are then solved for x and y using Gauss-Seidel iteration with Successive 
Over Relaxation (SOR). If the solution for x and y has not converged sufficiently then the 
interior point distribution just determined is used to again calculate P and Q on the 
boundary. The process is repeated until the smallest change between iterations for all values 

of x and y is below some threshold, e.g., 1.0 x 10 -4 . 

Figure 4. 12 illustrates the results of the modified control functions, Eqs. (4. 17) and 
(4.18), for two different cases of top and bottom boundary point distributions. The first 
case has the same point distribution on the top and bottom boundaries. The second has a 
top boundary point distribution that is clustered closer around the nozzle throat than the 
bottom boundary point distribution The vertical grid lines in Fig. 4.12, lines of constant 
can be seen to be orthogonal along the top and bottom boundaries. 

Figure 4.13 contains a magnified view of the upper boundary at positions upstream 
of the nozzle throat and near the nozzle exit. Comparison between Figures 4.9 and 4.13 
demonstrate the improved orthogonality at the boundary upstream of the nozzle throat using 
the modified control functions. However, at the nozzle exit, along the top boundary, it can 
be seen that the grid points in the interior have clustered very closely to the wall. Because 
the spacing of the right vertical boundary has not been maintained adequately in the interior, 
the grid lines, lines of constant T|, are unacceptably skewed. The close clustering near the 
wall was caused by the additional control function [46]. 
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Figure 4.12 Modified control function grid distributions for an 
arcjet thruster geometry (a) same point 
distributions, (b) different point distributions, on 
top and bottom boundaries 
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(b2) 


Figure 4.13 Magnified view of modified control function grid 
distributions for an arcjet thruster geometry (d) 
upstream nozzle wall, (2) nozzle exit) with, (a) 
same point distributions, (b) different point 
distributions, on top and bottom boundaries 
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To eliminate the skewed grid lines at the nozzle exit the spacing of the first grid 
point from the top and bottom boundaries was specified by Sorenson [48]. The spacing 
between the first grid point and the boundary on a vertical boundary line was multiplied by 
the ratio between the height of the vertical boundary and the spacing between the top and 
bottom boundaries on a line of constant £. The left boundary spacing was used for points 
upstream of the throat exit and the right boundary spacing was used for points downstream 
of the throat exit 

In addition, a limit was imposed on the change per iteration of the additional control 
functions shown in Eqs. (4.17) and (4.18). Sorenson’s [48] limiting function was 
modified to 

ZJujv = Z£urv + SIGNjmin [(0 p‘| Zeurv ” ZJurv I i Plim' max (|Zcurv| ,0.1 )] jZcurv " ^?urv) 

(4.19) 

where Z curv = — — — on the top and bottom boundaries and 

l r nl 2 

where = — — j — on the left and right boundaries. 

ZJjurv is the value of the curvature term from the previous, or nth iteration and Z^p, is the 
currently calculated value. Thus, the value of the curvature term to be used during the 
current iteration, ZJ&l, is damped by a combination of under-relaxing and limiting changes 
to a fraction of their currently calculated value. The SIGN function returns the magnitude of 
the first argument with the sign of the second argument Values of co p = 0.3 and pum= 0.01 
were used to obtain the results which follow. 

The grid generation algorithm was modified to allow the spacing of the first grid 
point away from the top and bottom boundaries to be fixed a priori. In addition, the 
limiting function, discussed above, was included to damp oscillations in the early stages of 
the solution. To aid in obtaining orthogonal grid in the interior of the region as well as at 
the boundary, only grids with a tighter top boundary point clustering were examined for 
use. 

Figure 4.14 contains the results for a top boundary point distribution more tightly 
clustered around the nozzle than the bottom boundary point distribution. The skewed grid 
was eliminated in the nozzle exit by including spacing control for the first grid point near 
the boundary. In addition, near orthogonality was maintained at the boundary throughout 
the region. 

Numerical grids used in the following computations are of three types. The most 
common grid is generated using the Poisson equation with the Thomas and Middlecoff 
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control functions and does not produce a grid orthogonal to the boundary. Comparison 
grids used to evaluate the effect of the numerical grid on the solution accuracy are generated 
using the Poisson equation with the Thompson control functions (with or without spacing 
control) and do produce a grid orthogonal to the boundary. Either the grid used will be 
displayed or the method used to generate the grid will be mentioned. 
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Figure 4.14 Modified control functions with spacing control 
grid distributions for an arcjet thruster geometry, 
different point distributions on top and bottom 
boundaries, (a) entire thruster, magnified views 
of (b) upstream nozzle wall and (c) nozzle exit 
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Chapters 


Results and Discussion 


5.1 Introduction 

In this chapter the algorithm used to solve the Thin Layer Navier-Stokes equations 
in transformed coordinates will be validated by comparing numerical predictions with 
experimental data. Available experimental results include Mach contours, wall static 
pressure ratios, centerline Mach number distributions, and integral parameters. 

Predicted results will illustrate the effects of geometry, viscosity, and swirl on the 
nozzle performance.The important geometries that must be accurately modeled include large 
area ratio nozzles (both inlet-to-throat and exit-to-throat area ratios) and nozzles with 
centerbodies. Combining these geometric effects with viscosity and swirl will enable 
swirling flow in arejet thrusters to be examined. 

Six nozzle geometries are examined. Hie first two nozzles were originally examined 
by Dutton [40, 41]. 1. Dutton's initial nozzle will illustrate the effects of viscosity and inlet 
swirl velocity profiles (Fig. 5.1) on nozzle performance parameters and Mach contours. 2. 
Dutton's experimental nozzle will demonstrate the effect of viscosity and inlet total 
pressure and swirl velocity profiles on static wall pressures, Mach contours, and nozzle 
swirl velocity profiles. 3. Back & Cuffel's nozzle [49-51] clearly illustrates a reflected 
oblique shock caused by the wall geometry. 4. An annular nozzle will show the effect of a 
viscous boundary layer and highly swirled flow in a nozzle containing a centerbody. 5. A 
100:1 area ratio nozzle containing a thick boundary layer will then be examined. 6. A 
converging-diverging nozzle containing a centerbody will be used to illustrate viscous, 
swirling flow through an arejet thruster geometry. 

5.2 Previous Numerical Research 

Early work in swirling flow was limited by the analysis approach chosen , quasi- 
one-dimensional [52] for example, or the use of a single level of swirl [53]. Later work 
examining swirling flow in annular nozzles [54, 55] showed the presence of swirl 
decreases the discharge coefficient, thrust, and vacuum specific impulse of a nozzle. 

Dutton's [40, 41] work in inviscid swirling flow used an explicit finite difference 
technique. Converging, converging-diverging, and annular nozzles were examined. The 
effects of three different inlet swirl velocity profiles on integral parameters were 


78 


determined. The results showed the degradation in performance due to swirl. A series of 
experiments were also conducted to compare measured and numerically predicted wall 
static pressures in a second converging-diverging nozzle. The numerical results showed 
excellent agreement with the experiments except in a region just downstream of the throat. 
As will be demonstrated shortly (see Section 5.4.2 entitled Dutton's experimental nozzle) 
the discrepancy was caused by the absence of viscosity. 

In all of the numerical investigations mentioned above, the flow was assumed to be 
inviscid. Chang [39, 42] was the first to examine viscous swirling flow in the three nozzle 
geometries used by Dutton. An implicit ADI technique was used in transonic regions. To 
eliminate numerical instabilities in the central-differenced ADI procedure in regions where 
the flowfield is predominantly supersonic, the procedure had to be modified. This 
modification was accomplished by retaining the central differencing in the cross- stream 
direction and flux-vector splitting in the stream wise direction. 

In a low Reynolds number flow, such as those found in arcjets [7], viscous effects 
cannot be neglected because of the large growth of the boundary layer in the nozzle exit 
[56]. A robust and efficient numerical technique must be used to solve the resulting 
governing equations. The present study was conducted in order to develop an efficient 
method of assessing the effects of viscous swirling flow in these nozzles. 


5.3 Nozzle Performance Parameters 

Several integral parameters will be used in assessing the effect of swirl on nozzle 
flows. The discharge coefficient, Cd, defined as 


C d h>- = 
m id 



(^t-^t)(p* u *)id 


(5.1) 


is a measure of the decrease in mass flow rate due to flow geometry and swirl. The vacuum 
stream thrust efficiency, Tivs. is defined as 


Tlvs 


Tid 


J prwe 

(P + pu 2 )r dr 

fee 

( r we " r ?e)(Pe Pc^ejid 


(5.2) 


the specific impulse efficiency, t]si, is a combination of the above two parameters and is 
defined 


Tlsi = 


(T/m) _ T| vs 

(T/mJtd C D 


(5.3) 
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%v s . 


Finally the nozzle inlet swirl number, Si, is defined as 

r 

I puwr 2 dr 

Si = — (5.4) 

f T wi 

r W i| pu 2 r dr 

* r ci 

and is the axial flux of angular momentum divided by the inlet wall radius times the axial 
flux of axial momentum. It is a direct measure of the level of swirl at the nozzle inlet. 

The integral parameters defined above arc non-dimensional by default.The values 
cited later for thrust and mass flow rate are also non-dimensional by virtue of being 
calculated with non-dimensional values of density, velocity, etc. 

5.4 Computations and Comparisons for a Variety of Nozzles 

A few words are needed regarding the conditions used to obtain the results in 
Sections 5.4.1 -5.4.6. As discussed in Section 3.3.2.2, the inlet boundary conditions 
specified at the subsonic inlet include the ratio of swirling velocity to axial velocity, <}>= 
w/u, stagnation pressure, P t , stagnation temperature, T t , and the ratio of radial velocity to 

axial velocity, 0= v/u. The radial velocity is always set equal to zero (0= v/u= 0.0) and 
unless otherwise stated the stagnation pressure and temperature are set equal to one. The 
remaining boundary condition for the swirl velocity component, w, is set equal to zero 
(<|)= w/u= 0.0) in the case of unswirled flow. The same inlet swirling velocity profiles used 
by Dutton, as shown in Fig. 5.1, are used for swirling flow cases here, except the swirl 
angle asymptotically approaches zero at the wall in all viscous calculations. 


Forced Vortex 



Figure 5.1 Inlet swirl velocity boundary condition profiles 
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A linear forced vortex profile has been assumed near the centerline for the constant 
angle and free vortex profiles to enforce the centerline boundary condition w= 0. The non- 
dimensional radius ratio of R/R W all = 0-2 was chosen as the matching point by Dutton 
because it was typical of geometries he examined. 

In all nozzles a constant specific heat of y 0 = 1.4 was assumed and Sutherland's 
formula [27] was used for the viscosity (constants were for air). 

5.4.1 Dutton's Initial Nozzle 

The 35*- 18.5* converging-diverging nozzle initially examined by Dutton [40, 41] 
has an exit-to-throat area ratio of 2.56. Two different numerical grids were used to examine 
the flowfield and are shown in Fig. 5.2. The grid used for the inviscid flow case had 45 
equally spaced grid points in the axial and 30 equally spaced grid points in the radial 
direction. From here on this will be written 45 x 30 grid points (axial x radial). The grid 
used for the viscous flow case had substantial clustering near the wall and contained 45 x 
40 grid points. 

Figures 5.3 and 5.4 illustrate the inviscid and viscous (Re= 10,000) Mach contours 
for both unswirled and swirling (30* constant angle profile) flows. The inviscid contours 
agree well with Dutton's results and the viscous contours agree well with Chang's results 
[42]. The main effect of swirl is to shift the centerline Mach contours upstream which in 
turn decreases the mass flow rate (see Table 5.1). A thin boundary layer can be seen along 
the wall in both viscous cases. 

The increased swirling centerline Mach contour is due to an increase in the axial 
velocity near the centerline in swirling flow. The increase in axial velocity causes a decrease 
in the pressure (Bernoulli effect). The numerical results show the static temperature is 
essentially constant across the nozzle inlet for both unswirled and swirling flows and only 
drops a few percent in swirling flow. This combination of falling pressure and constant 
temperature lead to a drop in the gas density. The drop in density is significant across the’ 
nozzle inlet and leads to the decreased mass flow rate. 

These effects are illustrated in Fig. 5.5 for the inviscid Mach contours shown in 
Fig. 5.3. The increase in axial velocity due to swirl is confined to a region less than 40% of 
the nozzle radius. However, the decrease in density due to swirl is much larger over a 
larger portion of the nozzle radius and hence leads to the decrease in mass flow rate shown 
in Table 5.1. 

It should also be noted that the converged numerical solutions, in both unswirled 
and swirling flows, exhibit an inlet axial velocity decreasing with increasing radius as 
shown in Fig. 5.5. Thus, the swirl velocity component for a constant angle inlet swirl 
boundary condition will not be constant in the outer portions of the nozzle inlet but will 
decrease proportional to the inlet axial velocity. 
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(a) 



Figure 5.2 Grid for Dutton's initial nozzle (DIN), (a) invisdd, 
45 x 30, (b) viscous, 45 x 40 
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Table 5.1 
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Nozzle performance summary for Dutton's 
initial nozzle 
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Figure 5.5 Comparison of inviscid unswirled and swirling 
(30° constant angle) nozzle inlet density and axial 
velocity profiles in DIN 


Two additional inlet swirl cases were also examined. The viscous (Re= 10,000) . 
Mach contours for a 52.5* forced vortex and a 40* free vortex are shown in Fig. 5.6. Again 
the Mach contours, in the inviscid regions of the flow, agree well with Dutton's results 
(Chang did not show any results for these cases).The resulting discharge coefficients, mass 
flow rates, etc. are included in Table 5.1. Examination of Table 5.1 reveals viscous effects 
and swirl decrease the mass flow rate and hence the thrust in a converging-diverging 
nozzle. The Cd and T| V s also are smaller for the viscous and swirling flow cases while rjsi 
remains constant. Dutton observed the same trends for inviscid, swirling flow. 

The swirl velocity profiles at the inlet, throat, and exit planes for the constant angle, 
forced vortex, and free vortex inlet swirl flows are shown in Figs. 5.7-5.8. The maximum 
throat swirl velocity for each case was used to normalize the swirl velocity. The form of the 
inlet swirl profile is clearly maintained throughout the nozzles. The decay of the swirl 
velocity in the boundary layer is also evident near the wall in the viscous cases. 
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A constant angle inlet swirl profile produces the largest amount of swirl at the inlet 
for a given maximum swirl angle, i.e., <J> max , (see Figure 5.9). Consequently, the largest 
reduction in mass flow rate, which occurs for the highest amount of swirl, occurs for a 
constant angle inlet profile. Figure 5.10 illustrates the reduction in mass flow rate for a 
given (Jimax, by non-dimensionalizing all mass flow rates with the unswirled, viscous mass 
flow rate from Table 5.1. The free and forced vortex inlet swirl velocity profiles would 
require a much larger <t> max to have as noticeable an effect as the constant angle profile. At a 
4> m ax of 30 degrees the constant angle inlet profile mass flow rate decreases 4% compared 
to the unswirled case. The free and forced vortex profiles cause approximately a 1% 
decrease in mass flow rate at a 4> max of 30 degrees. Therefore, when other nozzles are 
examined later in this study, only the effects of constant angle inlet swirl profiles will be 
illustrated. 

Figure 5.11 compares Co and ti vs using Dutton’s inviscid results and the present 
viscous results. At a Reynolds number of 10,000 Co decreases 1% and tj vs decreases 
approximately 2% from the inviscid result. Dutton observed a constant value of T|si= 
0.971. The present results also showed t|si to be constant with values of 0.957 and 0.955 
for inviscid and viscous (Re= 10,000) flows, respectively. 

Figure 5.12 illustrates the effect of viscosity on Co and Ti vs for unswirled and 
swirling flow. As the Reynolds number increases and viscosity decreases Cq and T| vs 
approach their inviscid limits. 

Convergence histories are displayed in Fig. 5.13 which were obtained as a function 
of viscosity and the inlet swirl condition using the inviscid and viscous grids shown in Fig. 
5.2. The norm of the residual is defined as £ii,i.j(5Q/Qj where ii represents the five 
conserved quantities shown in Eq. (2.49). Grid density slows the convergence rate the 
most. Including viscous effects in the computations slows the convergence rate by a 
smaller amount. Swirl has little or no effect on the convergence rate. 

Mass conservation through the nozzle, corresponding to the 300th iteration of the 
convergence history cases just discussed, is shown in Fig. 5.14. Mass conservation is 
almost identical in all six cases and has a maximum error of approximately ± 0.6%. In all 
of the nozzles examined in this study mass flow, being an integral parameter, typically 
reached its steady-state value by the time the residual had dropped approximately four 
orders of magnitude. Contours of nozzle properties are virtually unchanged after the four 
order drop in residual has occurred with changes in properties occurring in the third or 
fourth significant figure. Thus, these mass conservation trends are typical for the nozzles 
examined in this study. 
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Figure 5.7 Swirl velocity profiles at the inlet, throat, and exit 
for a 30’ constant angle inlet profile in DIN, (a) 
inviscid, (b) viscous (Re= 10,000) flow 
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Figure 5.8 
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Figure 5.9 Dependence of inlet swirl number in DIN on inlet 
swirl profile and maximum inlet swirl angle in 
viscous flow (Re= 10,000) 



(de *> 

Figure 5.10 Non-dimensionalized mass flow rate in DIN for a 
given inlet swirl profile and maximum inlet swirl 
angle in viscous flow (Re= 10,000) 
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Figure 5.11 
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Figure 5.13 Convergence in DIN as a function of viscosity, 
grid size, and swirl (Re= 10,000 in viscous flow) 
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Figure 5.14 Mass conservation in DIN (axial mass flow rate 
divided by inlet mass flow rate, Re= 10,000 if 
viscous flow) 
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5.4.2 Dutton’s experimental nozzle 

Dutton [40,41] also examined a second converging-diverging nozzle both 
experimentally and numerically. The 35*- 18.5* nozzle has an exit-to-throat area ratio of 
2.56. A 50 x 40 grid, see Fig. 5.15a, clustered near the wall was used in the present study 
to obtain the following viscous results. Dutton's measurements included inlet stagnation 
pressures, swirl velocity profiles as well as static wall pressures. The measured distribution 
of stagnation pressure (see Fig. 5.15b) and closest fit to the forced vortex, constant angle, 
or free vortex swirl profiles described previously were used as inlet boundary conditions. 
Mass flow errors in the present study were less than 1% and the total residual dropped four 
orders of magnitude in 250 steps. 

Shown in Fig. 5.16a are the Mach contours for unswirled viscous flow. Figure 
5.16b compares the experimental static wall pressure data for unswirled flow with the 
predicted inviscid and viscous (Re= 10,000) results. The present inviscid results are 
essentially the same as Dutton's inviscid results, which show a slight dip in the predicted 
static wall pressure downstream of the throat. However, the beginning of the viscous 
boundary layer just downstream of the throat causes the viscous static wall pressure to 
almost exactly follow the experimental data. 

Figures 5.17 and 5.18 contain the Mach contours and static wall pressure 
distributions for two viscous (Re= 10,000) swirling flow cases. The inlet swirl velocity 
profile for Fig. 5.17 which most closely matched Dutton's experimental inlet swirl data 
was a 40* forced vortex profile. A 50* constant angle profile was used in Fig. 5.18 to 
match Dutton’s second set of inlet swirl data. The swirl velocity profiles at the inlet, throat, 
and exit planes for the forced vortex and constant angle inlet profiles are shown in Fig. 
5.19. The maximum throat swirl velocity for each case was used to normalize the swirl 
velocity. The form of the inlet swirl profile is maintained throughout the nozzle and the 
decay of the velocity is clear in the boundary layer. 

In both swirling flow cases the inviscid results may be seen to underpredict the 
static wall pressure distributions downstream of the throat while the viscous result shows 
excellent agreement with the data. It should be noted that downstream of an axial location 
of approximately z= -0.5, the measured static wall pressure distributions for all three flow 
cases are essentially identical. Thus, the measured inlet stagnation pressure distribution has 
a large effect only in the inlet region where the swirl velocity is highest 

The values of Cq, T| V s. and T]si shown in Table 5.2 are 1-2% lower than Dutton’s 
reported values, exhibiting the degradation of performance due to viscosity. 

Comparing Tables 5.1 and 5.2 one sees the values of Sj for the forced vortex and 
constant angle cases are of the same order of magnitude for Dutton's initial and 
experimental nozzles. However, the Mach contours (Figs.5.6a, 5.17a and Figs. 5.4b, 
5.18a, respectively) are much different in the inlet portions of the nozzle. The nozzles have 


95 
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wall 


Figure 5.15 (a)ViscousgridinDutton'sexperimentalnozzle(DEN) 
(50 x 40) and (b) inlet stagnation pressure profiles for 
Dutton's experimental nozzle 
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O Dutton experimental data 
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Inviscid result 
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Figure 5.16 Viscous results for a 40* forced vortex inlet swirl 
profile in DEN, (a) Mach contours, (b) comparison of 
experimental and numerical static wall pressure dis- 
tributions (Re= 10,000) 
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O Dutton experimental data 
Viscous result, Re= 10,000 



Figure 5.17 Viscous results for a 50* constant angle inlet swirl 
profile in DEN, (a) Mach contours, (b) comparison of 
experimental and numerical static wall pressure dis- 
tributions (Re= 10,000) 
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Figure 5.18 Viscous results for a 40* forced vortex inlet swirl 
profile in DEN, (a) Mach contours, (b) comparison of 
experimental and numerical static wall pressure dis- 
tributions (Re= 10,000) 
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Figure 5.19 Swirl velocity profiles at the inlet, throat, and exit 
in DEN for a (a) 40* forced vortex, (b) 50° 
constant angle inlet profile in viscous flow (Re= 
10 , 000 ) 
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Table 5.2 


Nozzle performance summary for Dutton's 
experimental nozzle 


Figure 

Re 

Grid 

Swirl 

C D 


^SI 

St 

Thrust 


5.16 

10,000 

50x40 

— 

0.966 

0.918 

0.950 

0.0 

4.32 

2.40 

5.17 

10,000 

50x40 

m 

0.957 

0.910 

0.951 

0.355 

4.44 

2.47 

5.18 

10,000 

50x40 

Hi 

0.935 

0.888 

0.950 

0.766 

4.48 

2.49 







































the same exit-to-throat area ratios of 2.56 but vastly different inlet-to-throat area ratios: 
initial nozzle= 1.5 and experimental nozzle= 4.0. In a large inlet-to-throat area ratio nozzle 
the inlet axial velocity will be low, approaching stagnation conditions, making it easy to 
achieve a high ratio of swirl to axial velocity, w/u, and hence a large inlet swirl value, Si. 
While the ratio may be high the absolute magnitude of w is small causing little effect on the 
centerline axial velocity distribution or Mach contours in the inlet region. 

Comparing the Mach contours between the two nozzle geometries downstream of 
the throat reveals few differences. Since the throat area is much smaller than the inlet area 
the throat axial velocity is very high as the flow accelerates in the converging portion of the 
nozzle. Due to the exit boundary conditions the flow continues to accelerate in the diverging 
portion of the nozzle where the axial velocity component dwarfs the swirling velocity 
component. Due to conservation of angular momentum, the swirl velocity component is 
highest at the throat (smallest area) and lower as the area increases (see Figs. 5.7 and 
5.19). Thus, the swirl velocity in the diverging portion of the nozzle has a much lower 
magnitude than the axial velocity and consequently has a small effect of the Mach contours. 

Also worth noting in Table 5.2 is the increase in mass flow rate in the two swirled 
cases over the unswirled case. Typically swirl decreases the mass flow rate. However, the 
measured inlet stagnation pressure profiles. Fig. 5.15b, were used as boundary conditions 
for the unswirled and swirled cases. The swirled inlet stagnation pressures were higher 
than the unswirled case and more than offset the decrease in mass flow rate due to swirl for 
an overall increase in mass flow rate. 
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5.4.3 Back and Cuffel nozzle 


Back and Cuffel [49-51] performed a detailed series of measurements on a variety 
of converging-diverging nozzles. In one of the nozzles, the exit length was of sufficient 
length for an oblique shock to coalesce and reflect from the nozzle centerline. A 60 x 40 
grid (Fig. 5.20a) clustered near the wall was used to obtain the following viscous results 
for this nozzle (hereafter referred to as the B&C nozzle). The nozzle exit-to- throat and inlet- 
to-throat area ratios are 9.75 and 6.55, respectively. 

Figure 5.20b shows the Mach contours for a viscous flow with a Reynolds number 
of 10,000. These results were obtained for the nozzle examined experimentally by B&C. A 
growing boundary layer is visible along the wall in the diverging portion of the nozzle. 
Also visible is an oblique shock which forms downstream of the throat. 

The oblique shock is due to a discontinuity in the wall curvature. The wall geometry 
in this region is represented by a circular arc throat and a conical divergent wall. While the 
slope of the wall is continuous the wall curvature is discontinuous, leading to a 
compression of the flow downstream of this point. The beginning of the formation of an 
oblique shock can also be seen in the nozzles examined in Sections 5.4.1- 5.4.2. However, 
the divergent portion of the nozzles are too short for the shock to reflect from the centerline. 

Figure 5.20c contains the predicted and measured centerline Mach number 
distributions. The intersection of the shock at the centerline is predicted well and can be 
better resolved if more grid is used in the axial direction. 

Mach contours and centerline Mach number distribution for the present work are in 
good agreement with the inviscid results shown in Loth et al. [57] and Serra [58]. Loth 
used a finite element method with an adaptive grid scheme employing 4,990 elements to 
improve the shock capturing ability, while Serra employed a Lax-Wendroff procedure. 

Figure 5.21a illustrates the effect of a 60* constant angle inlet swirl profile on the 
viscous Mach contours. Identical contour levels (starting at M= 0.1 with 0.1 intervals) are 
shown in Figs. 5.20b and 5.21a. Again, as seen previously, swirl causes an upstream 
shifting of the Mach contours. It was found that swirl (Sj= 0.568) causes Co to drop 3.2% 
from 0.967 to 0.936, while the mass flow rate dropped 2.4% from 2.45 to 2.39. Thrust is 
similarly decreased by 3.3% from 4.80 to 4.64. As in the previously examined geometries 
T|si is again constant with a value of 0.944. 

Comparing the axial distribution of Mach number between the unswirled and 
swirled cases, Fig. 5.21b, shows a significant difference only in the inlet portion of the 
nozzle. As discussed previously, the large inlet-to-throat area ratio results in a large swirl 
number but relatively low absolute value of swirl velocity magnitude. Conservation of 
angular momentum results in a low swirl velocity in the divergent portion of the nozzle in 
comparison to the rapidly accelerating axial velocity. Hence, little change is observed in the 
Mach number distribution downstream of the throat. 
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Figure 5.20 Viscous results for Back & Cuffel (B&C) nozzle in 
unswirled flow, (a) 60 x 40 grid, (b) Mach contours, (c) 
comparison of experimental and numerical static wall 
pressure distributions (Re= 10,0002 
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Figure 5.21 Viscous results for Back & Cuffel nozzle, (a) Swirling 
Mach contours, 60* constant angle, S,= 0.568, (b) 
swirled and unswirled centerline Mach number distri* 
butions (Re= 10,000) 
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The swirl velocity profiles at the inlet, throat, and exit planes for the 60* constant 
angle inlet swirl profile are shown in Fig. 5.22. The classic constant angle profile shape is 
clearly visible in the inlet and throat profiles and to a lesser degree in the exit profile. The 
large inlet-to- throat area ratio combined with conservation of angular momentum causes the 
magnitude of the throat profile to be the largest 

The residual for the unswirled case dropped six orders of magnitude in 500 
iterations while the residual in the swirled case only dropped four orders of magnitude in 
750 iterations. Clearly swirl has a significant effect on the convergence rate. This is 
contrary to the finding for Dutton's initial geometry. The reason for the decrease in 
convergence rate is two-fold. 

First, the inlet-to-throat area ratio for the B&C nozzle is 6.55 while it was only 1 .5 
for Dutton's initial geometry. A quasi one-dimensional approximation is used as the initial 
flowfield profile. Clearly, small inlet-to-throat area ratio nozzles have converged solutions 
which deviate little from this initial guess because the flow is almost one-dimensional. As 
the inlet-to-throat area ratio increases the final converged solution for the flowfield is highly 
two-dimensional. Thus, a highly two-dimensional flowfield requires more' iterations to 
reach a converged solution than does a flowfield approaching a one-dimensional 
approximation. 

The second reason the convergence rate is decreased is because the large increase in 
inlet area occurs in a region of subsonic flow. Information can propagate both upstream 
and downstream in a subsonic flowfield. Thus, the inlet region of the nozzle requires many 
iterations to reach a steady-state solution. The supersonic portion of the nozzle can only 
reach convergence after the subsonic portion converges because the flow of information in 
a supersonic flowfield is only in the downstream direction. 

Checking conservation of mass revealed a maximum error of 1.9% just upstream of 
the throat for the 60 x 40 grid (for both unswirled and swirled flow cases). A grid 
refinement study was undertaken to examine the effects of the number of grid points on the 
converged swirling solution. Running successively coarser grids of 45 x 30 and 30 x 20 
increased the maximum mass flow error to 3.1% and 4.6%, respectively. The intersection 
of the oblique shock on the centerline moved downstream with decreasing grid density. 
Finally, a 90 x 40 grid was run resulting in a maximum mass flow rate error of less than 
1% and a centerline/oblique shock location very slightly upstream of the 60 x 40 location. 
Thus, the solution approaches an asymptotic limit as the number of grid points increase. 

A 60 x 40 grid with the cross-stream mesh lines orthogonal to the nozzle wall (Fig. 
5.23) was also run. Integral parameters varied less than 2% from the results obtained with 
the conventional grid shown in Fig. 5.20a and the Mach contours are essentially unchanged 
from the results shown in Fig. 5.21b. The strong conservation form of the governing 
equations adequately represents the flowfield irregardless of the numerical mesh providing 
a sufficient number of grid points is used. 
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Figure 5.22 Swirl velocity profiles at the inlet, throat, and exit in 
B&C nozzle for a 60° constant angle inlet profile in 
viscous flow (Re= 10,000) 



Figu re 5.23 Viscous grid in B & C nozzle, 60 x 40, with cross-stream 

mesh lines orthogonal to nozzle boundaries 
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5.4.4 Annular Nozzle 


An annular nozzle, also examined numerically by Dutton [40, 41] and Chang [39, 
42], was examined in order to investigate nozzles with centerbodies. The upper surface is a 
constant radius wall. The lower surface has a cylindrical inlet with a circular arc transition 
to a circular arc throat and 10* conical wall. A 50 x 40 grid (Fig. 5.24a) clustered near the 
upper and lower surfaces was used to compare to Dutton's inviscid results. The unswirled 
viscous Mach contours (Re= 10,000) are shown in Fig. 5.24b. The inviscid portion of the 
unswirled flowfield agrees well with Dutton's inviscid result, but significant boundary 
layers are evident on both surfaces. The present viscous, unswirled results agree almost 
identically with Chang’s work. 

A 70* constant angle inlet swirl profile (maximum swirl velocity, w, is =2.75u, 
thus the inlet flow is primarily tangential!) was used to obtain the swirling results shown in 
Fig. 5.25a. As observed previously swirl causes the Mach contours to shift upstream.The 
resulting inlet Mach contours are much different than Dutton’s inviscid result because of the 
developing boundary layer on the lower surface. However the viscous, swirling Mach 
contours are again in almost identical agreement with Chang's results. 

Swirl (Sj = 1.85) causes Co and the mass flow rate to drop 9.9% from 0.923 to 
0.832, and from 0.973 to 0.875, respectively. Thrust is similarly decreased by 11.0% 
from 1.63 to 1.45. In addition, rjsi is not constant as in the earlier geometries but decreases 
1.5% from 0.989 to 0.974 in swirled flow. These decreases in nozzle performance were 
the largest calculated in this study. 

It should also be noted that an extremely large swirl angle was needed as an inlet 
boundary condition because of the large inlet-to-throat area ratio of 3.06 (the exit-to-throat 
area ratio is 1.37). The resulting axial velocity component is low at the inlet because of the 
large area. Inlet swirl angles less than 30* produced values of Si less than *1.0 and the 
resulting flowfield is hardly altered from the unswirled case. Hence, a very large inlet swirl 
angle is required to produce a noticeable effect on the Mach contours since a large value of 
w/u is easy to achieve given the small inlet value of u. 

The swirl velocity profiles at the inlet, throat, and exit are shown in Fig. 5.25b. 
These profiles follow a much different trend than observed earlier in the converging- 
diverging nozzles without centerbodies. In the inviscid core portion of the flowfield, i.e., 
outside the viscous wall layers, the profiles mirror Dutton's results. A constant angle 
profile is only apparent in the inlet profile, while the throat and exit profiles are essentially 
flat outside of the boundary layer. Apparently the reduction in area coupled with a relatively 
small radius change across the throat and exit planes of the nozzle causes the angular 
momentum to be essentially constant across the cross-section of the nozzle. Angular 
momentum appears to be preserved since the relative magnitudes of the swirl velocity are 
inversely proportional to the centerbody radius at the inlet (smallest), throat (largest), and 
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Figure 5.24 Viscous results for an annular nozzle (a) 50 x 40 grid, 
(b) unswirled Mach contours (Re= 10,000) 
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w/max(w |nlet ) 


Figure 5.25 Viscous results for an annular nozzle, (a) swirling 
Mach contours, 70* constant angle (Si= 1.85), (b) swirl 
velocity profiles at the inlet, throat, and exit (Re= 
10 , 000 ) 
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exit (intermediate). These results could not be contrasted to Chang as he did not present any 
results for the swirl velocity components. 

An exact value of the swirl or axial velocity cannot be specified at the inlet, only the 
ratio of swirl/axial. As stated earlier, in viscous, swirling flow cases the inlet swirl angle is 
modified to asymptotically approach zero at the wall. This modification to the constant 
angle inlet swirl profile was enforced at the upper surface of the annular nozzle (Fig. 5.25). 
It was realized that this modification should also be enforced at the surface of a centerbody 
also. It is clear from Fig. 5.25b that the axial velocity does not approach zero as the 
centerbody surface is approached since the normalized swirl velocity component is * 0.9. 

Figure 5.26 (Sj= 0.987) illustrates the drastic modification to the flowfield when 
the inlet swirl angle is modified to asymptotically approach zero at both the upper and lower 
solid surfaces. The Mach contours are not shifted as far upstream as in Fig. 5.25a and the 
growth of the boundary layer occurs quickly but with a smoother transition above the 
centerbody. The swirl velocity profiles (Fig. 5.26b) show the same general trends in the 
inviscid region as for the prior case (Fig. 5.25b) but the boundary layer development is 
now much clearer on both surfaces as the swirl velocity approaches zero at a solid surface. 
Apparently Chang only used an asymptotic swirl velocity profile at the upper solid surface 
since Fig. 5.25a almost identically matches his viscous swirling flow results. 

Values of Cq, mass flow rate, thrust, and “Hsi are 0.882, 0.928, 1.54, and 0.974, 
respectively, for the asymptotic inlet swirl case at both upper and lower solid surface. Table 
5.3 contains a summary for the annular nozzle cases discussed in this section. 

Examination of mass conservation revealed there was a maximum error of 0.3% 
based on the inlet mass flow rate in the three cases discussed above. The residual dropped 
four orders of magnitude in 500 iterations for the unswirled case and only two orders of 
magnitude in 500 iterations and three orders of magnitude in 750 iterations for the swirled 
cases. However, there was less than 1% change in the swirl case integral parameters 
between 500 and 750 iterations and the Mach contours showed only slight changes. For 
qualitative purposes the solution had converged in 500 iterations, but, if detailed flowfield 
information was desired, the swirled cases would need to run longer until the residual had 
dropped at least four orders of magnitude. 
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Figure 5.26 Viscous results for an annular nozzle with a modified 
asymptotic inlet swirl profile, (a) swirling Mach con- 
tours, 70* constant angle (Si= 0.987), (b) swirl velocity 
profiles at the inlet, throat, and exit (Re= 10,000) 
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Table 53 


Nozzle performance summary for the annular nozzle 


Figure 

Re 
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0.875 
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5.4.5 100:1 Area Ratio Nozzle 


The largest area ratio nozzle examined to present using the methodologies 
developed in this work is a 100:1 exit-to-throat nozzle used in experimental resistojet 
studies at NASA Lewis Research Center [59]. The experimental inlet stagnation pressure 
and temperature were 6,400 Pa and 699 K, respectively, with a throat Reynolds number of 
850. A 60 x 45 grid with clustering near the wall was used to obtain the numerical results 
for the 45*-20* nozzle shown in Fig. 5.27a. Also included in the figure, at the same scale, 
is the numerical grid from the experimental converging-diverging nozzle (50 x 40 grid) 
examined earlier. 

Figure 5.27b, c contain the viscous Mach contours for unswirled and swirling 
viscous (Re= 850) flows, respectively. In both cases the boundary layer occupies a large 
portion of the diverging section of the nozzle. It should be noted that the exit plane 
centerline Mach number is approximately 5.0 in each flow. Figure 5.28 contains the Mach 
contours for an unswirled viscous flow at Re= 5,000. Except for a thinner boundary layer 
and a slight shift downstream of the oblique shock/centerline intersection point the 
flowfield looks much the same as the un swirled, Re= 850 case. 

Figure 5.29 compares measured and predicted exit plane pitot pressures. The 
predicted result at Re= 850 is in good agreement with measurements near the centerline but 
differs with the experimental data as the radius increases. The pitot pressure profile for the 
swirling flow case at a Re= 850 was essentially the same as the unswirled result. Several 
reasons for the disparity between the predicted and experimental results may be offered. 

One possible reason is inadequate grid density. As may be concluded from 
examination of Fig. 5.27a, a 60 x 45 grid is probably the minimum size needed to obtain 
adequate resolution of the flowfield for this large area ratio nozzle. But as seen in the 
oblique shock nozzle case discussed earlier, inadequate grid density leads to mass 
conservation errors. The residual for the unswirled case dropped 4 orders of magnitude in 
1,300 iterations, indicating the flowfield would change little from the current values, but 
had a maximum mass flow rate error (based on the inlet mass flow rate value) of 16% in 
the subsonic region and a constant 10% error downstream of the throat. Thus a converged 
solution with poor mass conservation points to inadequate grid density (see the B&C 
nozzle discussion). 

The probability that inadequate grid density in the subsonic portion of the nozzle 
leads to mass conservation errors is large when one compares the inlet-to-throat ratios of 
the two nozzles in Fig. 5.27a. The inlet area ratios of the Back & Cuffel/oblique shock 
nozzle and the 100:1 area ratio nozzles are 9.75 and 48.4, respectively. The nozzles have 
roughly the same number of grid points in the radial direction (40 and 45 respectively). The 
B&C nozzle required 60 grid points in the axial direction to reduce the mass conservation 
error to less than 1%. Approximately 20 grid points were in the subsonic portion of the 
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Back & Cuffel 




(c) 



Figure 5.27 Viscous results for a 100: 1 area ratio nozzle (a) 60 x 45 

grid, (b) unswirled Mach contours, (c) swirling Mach 
contours, 30 s constant angle (Re= 850) 
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Figure 5.28 Unswirled Mach contours for a 100:1 area ratio nozzle 

(Re= 5,000) 


« Experimental data 



0.0 4.0 8.0 12.0 16.0 

Radius (mm) 


Figure 5.29 Comparison of measured and predicted exit plane 
pitot pressure distribution in 100:1 area ratio nozzle 
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nozzle. The 100:1 area ratio nozzle employs approximately 25 grid points in the axial 
subsonic portion of the nozzle. Thus roughly the same number of grid points were used to 
represent areas a factor of five different 

Another source for error is an insufficient number of iterations. The residual for the 
swirled flowfield had only dropped one order of magnitude in 1,300 iterations. Checking 
the subsonic flowfield revealed the inlet plane had converged at 1,300 iterations but the 
remainder of the subsonic portion of the nozzle was still evolving. The mass flow rate was 
again constant in the diverging portion of the nozzle, indicating the supersonic portion of 
the flowfield was converged and changing as the upstream conditions evolved. 
Examination of the swirl profiles showed the swirl velocity at the throat and exit planes was 
still developing. 

Thus a combination of inadequate grid density and an insufficient number of 
iterations undoubtedly contributed to the poor exit plane pitot pressure prediction. 

Another possibility of error related to inadequate grid density lies with the oblique 
shock representation. The hump in the numerical pitot pressure prediction coincides with 
the position of the oblique shock reflecting downstream from the centerline to the exit 
plane. The reflected oblique shock could be inadequately represented because of the 
relatively sparse grid in the diverging portion of the nozzle. 

A final source for the error in the pressure discrepancy lies with the nozzle 
geometry. The throat radius of this nozzle was 0.0625" (0.0159 mm). However because of 
the machining process used and the extremely small dimensions of the nozzle the actual 
throat geometry, i.e. throat radius, radius of curvature of the nozzle walls, was not 
available at the time this study was performed. 

When an adequate model of the nozzle throat geometry is determined grid 
refinement studies will be conducted to obtain converged solutions with mass conservation 
errors of less than 1%. Then a comparison between predicted and experimental results can 
be made and the effects of swirl on the flowfield can be determined. 


117 



5.4.6 Arcjet thruster nozzle 

An arcjet thruster geometry combines all of the difficult geometrical features 
examined in the previous sections into one geometry. A typical arcjet geometry contains a 
centerbody, representing the cathode, that partially extends into a large area ratio 
converging-diverging nozzle. A geometry representative of arcjets [25] was formed by 
placing a centerbody with a 45' tip angle in the subsonic portion of a 35'- 18.5' 
converging-diverging nozzle. The C-D nozzle wall contour is the same geometry used in 
Dutton's experimental nozzle examined earlier. The inlet and exit to- throat area ratios are 

3.6 and 4.0, respectively. 

A 60 x 45 grid shown in Fig. 5.30a was clustered at both surfaces and was used to 
resolve the boundary layers on the centerbody and nozzle wall. Figure 5.30b shows the 
resulting Mach contours for viscous flow at a Re= 850. As could be expected from the 
previous results a thick boundary layer forms on the nozzle wall because of the very low 
Reynolds number. Viscosity obviously plays an extremely important role in flow. The flow 
rapidly accelerates downstream of the cathode tip along the nozzle centerline. However, the 
sonic line extends far downstream of the nozzle wall sonic line location. Comparing the 
location of the centerline sonic point with the result for Dutton's experimental nozzle shows 
the centerbody has pushed the sonic line further downstream. Although this is not a 
surprising result, it demonstrates the highly two-dimensional nature of the flow. 

A 30' constant angle inlet swirl profile modified to asymptotically approach zero at 
both surfaces (see Section 5.3. Annular nozzle for a discussion of the modification) was 
used to introduce swirl into the flowfield. The resulting viscous Mach contours are shown 
in Fig. 5.31a for a Reynolds number of 850. As observed previously the Mach contours 
are shifted upstream slightly as a consequence of the swirling velocity component. The 
resulting swirl number, Sj= 0.320, indicates a low level of swirl. Using a maximum swirl 
angle greater than 30* would have had a greater effect on the flowfield. 

Examining the integral parameters illustrates the small changes caused by swirl. The 
Cd decreased 1.2% from 0.896 to 0.885 while the mass flow rate decreased 1.3% from 
2.27 to 2.24. Thrust decreased 1.2% from 4.06 to 4.01 and T|si was constant at 0.946. 

The swirl velocity profiles at the inlet, throat, and exit are shown in Fig. 5.31b. The 
constant angle profile is visible in the inlet profile becomes distorted as the flow travels 
downstream from the inlet and is first compressed as it nears the cathode tip, expands over 
the cathode tip, and continues to expand in the diverging portion of the nozzle. 
Qualitatively, conservation of angular momentum is satisfied by noting the magnitude of 
the largest swirl velocity component occurs at the throat because the throat corresponds to 
the smallest area with the smallest mean radius. The exit- to- throat area ratio of 4.0 is 
approximately the same as the inlet- to- throat area ratio of 3.6. However, the mean radius of 
the inlet is larger than the mean radius of the exit plane. Thus the swirl magnitude is 
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Figure 5 JO Viscous results for an arcjet thruster nozzle (a) 60 x 45 

grid, (b) unswirled Mach contours (Re= 850) 
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Figure 5.31 Viscous results for an arcjet thruster nozzle, (a) swirl- 
ing Mach contours, 30* constant angle (Si= 0.320), (b) 
swirl velocity profiles at the inlet, throat, and exit (Re= 
850) 


120 


observed to be smaller at the inlet than at the exit plane. 

The maximum mass conservation error, based on the inlet value of mass flow rate, 
was less than 0.6% in both the unswirled and swirled cases. The residual for the unswirled 
case dropped five orders of magnitude in 1300 iterations, while the residual for the 
unswirled case only dropped three orders of magnitude in 1300 iterations. Thus, as 
observed earlier, swirl does affect the convergence rate in highly two-dimensional 
flowfields. 

Comparing the B&C and arcjet nozzles shows vastly different convergence rates for 
grids containing approximately the same number of grid points (« 60 x 40). In unswirled 
flow computations, the residual in the B&C nozzle case dropped six orders of magnitude 
while the residual in the arcjet case only dropped three orders of magnitude. The reason for 
the discrepancy is because the arcjet grid was clustered at both upper and lower surfaces 
while the grid used in the B&C nozzle was only clustered aldng the upper surface. Grid 
clustering causes a large decrease in the convergence rate. 

5.4 Computational requirements 


The data processing rate (DPR) for the viscous flow results using the implicit 
algorithm was 1.06 x 10*^ seconds/(grid point*iteration) on a Cray XMP 2/8 using default 
vectorization. A fully vectorized three-dimensional version employing MacCormack's 
algorithm reported [33] a DPR of 3.1 x 10' 4 on a Cray-2 single processor. Since the Cray 
2 is on the order of two times faster than a Cray XMP and little work has been done to 
vectorize the present code; it was felt that the DPR for this algorithm is acceptable. Recent 
use of Cray Research Inc. vectorization software indicates the DPR can be reduced by a 
factor of two with some minor modifications to the block-tridiagonal solver subroutines 
used in the code. 

As discussed in Section 3.3.1 for a grid described as 60 x 45, the computational 
domain actually contains 58 x 43 grid points for a total of 2,494 grid points. Thus, a 
solution for the arcjet nozzle requiring 1300 iterations would take »3,400 seconds of CPU 
time. 

The time step definition used by MacCormack [20] was modified for the non- 
dimensional form of the governing equations used in this study. The resulting definition for 
the time step based on a given CFL (Courant-Friedrichs-Lewy) number is 


CFL - 
At 


flu 1 d, + |v 1 d,)J + y ~ (dj + d£Jl : 
■ 2 ^ (dj 1 2d,d„ + dJ)j J 


p Re 
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The majority of computations in this study were run with an initial CFL value of 0.25 
incrementing to 25.0 by 0.25 increments/iteration. When the residual had dropped to one- 
tenth of the initial value, the CFL was incremented again by 0.5 increments/iteration until a 
maximum CFL value of 80.0 was reached and remained constant for the remainder of the 
computation. The swirled arcjet thruster geometry case was run with a maximum CFL of 
50.0. Higher maximum CFL numbers resulted in growing instabilities occurring at the 
nozzle centerline upstream of the exit plane. The instabilities were not seen in any of the 
other cases examined. 
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Chapter 6 


Summary & Recommendations 


6.1 Summary 

In this dissertation a method has been presented for modeling swirling flow through 
a constricted arcjet thruster geometry. The method allows the effects of swirl and viscosity 
to be determined in an axisymmetric nozzle (with or without a centerbody). The numerical 
method is stable for large time steps allowing steady-state solutions to be achieved in a 
small number of iterations. 

An order of magnitude analysis was used to examine all of the terms contained in 
the equations governing the fluid mechanics and electromagnetic fields in an arcjet thruster. 
Retaining the significant terms revealed that the remaining fluid mechanics and combined 
electromagnetic field equations are uncoupled. The only link between the two equation sets 
is through the electrical conductivity of the gaseous propellant, which is a function of the 
pressure and temperature. 

Experimental results have shown a swirling velocity component stabilizes the 
constricted arc. Thus, a circumferential or swirl velocity component as well as the axial and 
radial velocity components were needed in the analysis. Since a typical nozzle geometry is 
symmetric about the centerline, a two-dimensional axisymmetric formulation was utilized. 

An algorithm was developed to solve the governing fluid mechanics equations, the 
axisymmetric Thin-Layer Navier Stokes equations. Based on MacCormack's [20] 
algorithm, the technique is fully implicit, second-order accurate in space and flux-split By 
employing an implicit algorithm large time-steps could be used to reach the desired steady- 
state solution. Second-order spatial accuracy is required to obtain solutions valid 
throughout flowfields with complex geometries. Flux-splitting the inviscid fluxes 
accurately represents the physics of the flow by allowing information to travel in the 
appropriate direction. Thus, subsonic and supersonic regions are represented without the 
use of artificial damping. Since the method does not employ approximate factorization no 
additional terms are introduced into the finite-volume representation of the governing 
equations. 

An implicit boundary condition technique was successfully adapted for use at the 
inlet and exit nozzle planes. The method can be used to specify the appropriate number of 
conditions for subsonic or supersonic flows. The effects of three different inlet swirl 
velocity profiles were investigated by specifying different inlet plane ratios of swirl velocity 
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to axial velocity. Each inlet swirl profile produced a different effect on the flowfield and 
degraded the nozzle performance by varying degrees. 

A numerical grid generator was written to provide the mesh to represent the nozzles 
examined. The grid generator can produce mesh with the cross-stream grid lines non- 
orthogonal to the upper and lower surfaces, i.e., vertical, or with cross-stream grid lines 
orthogonal to the upper and lower surfaces. Since the governing equations were formulated 
in a strong conservation form, grid orthogonality should not affect the conservation of 
mass, momentum, and energy of the algorithm, but this could not be determined a priori. 
Grid orthogonality was found to have little effect on the converged steady-state solution. 

Since typical arcjets employ large area ratio converging-diverging nozzles 
containing centerbodies, a wide variety of nozzle geometries were needed to validate the 
algorithm. These included conventional converging-diverging nozzles with exit-to-throat 
area ratios as large as 100 (with inlet-to-throat area ratios as large as 48) and an annular 
nozzle with a full-length centerbody.Computed Mach number distributions, static wall 
pressures, and oblique shock structures were in excellent agreement with experimental data 
and previous numerical results. The degradation of nozzle performance due to swirl and 
viscosity was also clearly illustrated. Finally, a typical arcjet thruster geometry was 
examined. The effects of swirl on nozzle performance and of geometry on the swirl profiles 
throughout the nozzle were demonstrated. 

6.2 Conclusions 

The numerical algorithm provides an accurate method of calculating the viscous, 
swirling flowfields through converging-diverging nozzles. In both unswirled and swirling 
flows, comparisons with experimental data and other numerical calculations were in 
excellent agreement for quantities such as Mach contours, static pressure distributions and 
oblique shock structures including reflected oblique shocks. 

Viscosity was shown to play an important role in the wall pressure recovery 
downstream of the nozzle throat. Viscous effects were also shown to be important in large 
area ratio nozzles where the Re is on the order of 1,000. The boundary layer occupies a 
significant portion of the diverging portion of the nozzle in low Reynolds number/high area 
ratio nozzles. Viscous effects were also evident in nozzles containing centerbodies. The 
mass flow rate was also shown to decrease when viscous effects were included and 
compared with inviscid results. 

Swirl was shown to cause a decrease in the mass flow rate, and a resultant decrease 
in the thrust, through a rise in the centerline axial velocity distribution and a concomitant 
decrease in the density. This effect is observable through the marked upstream shifting of 
the Mach contours in highly swirled flows. The decreased mass flow rate can be offset by 
an increase in the throat radius or inlet stagnation pressure. Since the assumed inlet swirl 
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profile and maximum swirl angle greatly affects the resulting flowfield both need to be 
measured experimentally in order to obtain an accurate representation of the flowfield. 

The swirl velocity magnitude at a cross-section of a nozzle varies inversely to the 
area and mean radius of the cross-section (in comparison to the inlet area and inlet mean 
radius). As the flow accelerates through a nozzle throat into the diverging exit the ratio of 
swirl/axial velocity decrease due to the increase in area and axial velocity magnitude. Thus, 
the effects of swirl are most strongly felt in the subsonic portion of the nozzle. For instance 
very high inlet swirl levels would be needed to shift a centerline/oblique shock intersection 
point upstream in the diverging portion of a nozzle. 

Cross-sectional area magnitudes also strongly influence the convergence rate, 
especially in the subsonic regions of swirling flows. As the inlet-to-throat ratio increases 
the subsonic flowfield requires a longer time to converge to a highly two-dimensional 
flowfield from the one-dimensional approximation used as the initial guess. Since the initial 
guess assumes the swirling velocity component is zero the algorithm must determine the 
value throughout the flowfield. In highly swirled flow with a large inlet-to throat area ratio 
the change in the swirling velocity magnitude through the flowfield will be large. Because 
information propagates in all directions in the subsonic regime, convergence is rapidly 
slowed as the area ratio and level of swirl increases. 

Mass conservation is also affected by large inlet-to-throat area ratio nozzles. As the 
area ratio increases, a larger number of grid points is required to maintain mass flow rate 
errors less than 1%. With but one exception, the nozzles examined in this research had 
mass conservation errors, based on the inlet mass flow rate, of less than 0.6%. The 
exception was the 100:1 exit-to-throat area ratio nozzle which had an unconverged mass 
flow rate error of 16% in the subsonic inlet (inlet-to-throat area ratio of 48.4). More grid 
points and more iterations are required to achieve mass flow rate errors of less than 1% in 
the subsonic portions of large inlet-to-throat area ratio nozzles. 

A double penalty in convergence rate is paid to resolve the flowfield gradients in a 
large inlet-to-throat area ratio nozzles with swirl. More grid points are needed to resolve the 
large gradients but as the number of grid points increases the convergence rate decreases 
for any nozzle configuration. And as the inlet-to-throat area ratio increases the more the 
final converged swirling flowfield will differ from the initial assumed state of zero swirl, 
thus slowing the convergence rate. 

Convergence is not only affected by the number of grid points but by the location of 
the grid points. Clustering the grid points near a wall causes large decrease in the 
convergence rate for instance. Thus as the grid density increases the convergence rate 
decreases. Brief studies of grid orthogonality at the upper and lower boundaries indicated 
little change in the final solution but did decrease the increase the convergence rate slightly. 

The DPR of 1.06 x 10"^ seconds/(grid point*iteration) is adequate but leaves room 
for improvement. As more realistic arejet geometries are examined a faster processing rate 
becomes more important since the inlet- and exit-to-throat area ratios of arejets are large. As 
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discussed above the convergence rate decreases as the area ratio increases. 

6.3 Recommendations for Future Research 

The formulation presented in this dissertation showed the fluid mechanics and 
electromagnetic (Emag) field equations were decoupled, allowing the research to be 
conducted in two phases. The first phase has been completed by developing an algorithm to 
solve the governing fluid mechanics equations. The second phase would be to include the 
Emag field effects in an operational arcjet. This can be accomplished by developing an 
algorithm for the Emag field governing equation presented in Chapter 2. A solution for a 
running arcjet would be obtained by solving the fluid mechanics and Emag algorithms 
iteratively until a converged solution is reached. 

An alternative method for the second phase would be to reformulate the fluid 
mechanics algorithm to include the Emag field equation and solve the combined set of six 
governing equations simultaneously. 

At the present time a far reaching extension of the modeling efforts would be to not 
make an assumption of LTE. The resulting two temperature plasma model would include 
electron temperatures as well as ion and neutral particle temperatures. Thus, the model 
would be valid in regions near the cathode and anode surfaces where the single fluid 
plasma model of this research is not valid. 

Computational studies need to be conducted to resolve the mass conservation errors 
in the large inlet- and exit-to-throat area ratio nozzles. Obvious avenues of investigation 
including adding more grid points and/or running for more iterations. Partitioning the 
flowfield into sub and supersonic portions would be another approach. First, a converged 
solution should be obtained in the subsonic portion of the nozzle. The supersonic portion 
of the nozzle could then be solved using the subsonic portion as an inlet boundary 
condition. In this way computational time is not wasted on the supersonic portion of the 
nozzle when the subsonic portion, i.e., inlet condition, of the flowfield is still evolving. 

Another point needing resolution is the discrepancy in the exit plane pitot pressure 
distribution for the 100:1 area ratio nozzle. A more accurate nozzle geometry needs to be 
determined and combined with the computational studies above to provide a converged 
solution with mass conservation errors of less than 1%. 

Parametric studies could be conducted to determine the effects of nozzle shape on 
thruster performance. Alternatively, the fluid mechanics algorithm could become a portion 
of an algorithm to compute the optimum nozzle shape based on design parameters such as 
thrust, mass flow rate, and geometry. 

Other extensions to the algorithm would be to include an adaptive grid scheme to 
provide automatic resolution of high gradient regions or to include all of the viscous terms 
in the explicit portion of the algorithm, resulting in the full Navier-Stokes equations. 


Any large programming project will take twice as long 
as you estimate, even if you include that eventuality in 
your estimate. 

Babbage’s Law 

in honor of Charles Babbage (1792-1871 ), 
pioneer computer scientist 
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Appendix A 


Derivation of 

Quasi-neutrality Condition 


The justification for assuming a plasma is quasi-neutral arises from an examination 
of one of Maxwell's Equations: 

A 

<K-V — • 

V -E = p^/ej) (Al) 

where: A indicates a dimensional quantity 

p e is the net space charge density (coul/m 3 ) 
and Eq is the permitivity of vacuum (farad/m) 

=8.85 xlO'* 2 coul 2 sec 2 / (kg m 3 ) 

The following non-dimensional quantities are defined : 


V 


V 

f[h 


Pe = pe (e n e ) 


where e is the electron charge 

= 1.60207 xl0'19 coul / electron 

iig is the number of electrons / m 3 


and 

E = E (kg T / e r(h) 

kg is Boltzmann's Constant 

= 1.381 xlO' 23 (kg m 2 )/(sec 2 K) 

The quantity (kg T / e r^) is the ratio of the thermal energy to the electric charge of 
an electron and is a measure of the ability of an electron to move away from an ion. If the 
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distance an electron is separated form an ion is small then the plasma is essentially neutral 
in comparison to the macroscopic length scale of the flow, r^. The thermal energy (voltage 
potential) is not large enough for the electron to overcome its attraction to the ion and thus 
the flow is quasi-neutral over the macroscopic dimensions of the flow. However, the flow 
is not neutral on a length scale that has dimensions on the order of the electron-ion 
separation distance. 

Substituting the non-dimensional quantities into Eq. (Al) yields: 


V 



(A2) 


Recalling that the Debye shielding length, Xq, is the maximum distance over which 
an ionized gas may be non -neutral, on the average, depends on the number density of 
charge particles and their mean thermal speed and is defined as 


or 



£o ks T 

n e e 2 



(A3) 


Substituting Eq. (A3) into (A2) and rearranging yields 


V E V E 


/ItM 2 

r* 2 | 

(J 

i69 2 T / n e / 


(A4) 


For % = 2.0 xl0‘4 m, T = 10,000 K, and n e = 2x10^1 electrons/ m^ [26] we find: 


p c = ( 6.0x10 7 ) V E -0 
and 

X D = 1.54xl0" 7 m 


where Xq is the maximum distance over which an ionized gas may be considered non- 
neutral. Compared to the macroscopic length scale of the thruster, r t h, Xq is 10^ smaller. 
Consequently, p e , which is a measure of the n£i spatial charge of the plasma, is 
approximately zero and an assumption the plasma is quasi-neutral is justified. Thus, 

V E - 0 . 
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Appendix B 


Coordinate Transformation 
and 

Derivation of the 

Transformed Governing Equations 


Consider the transformation from the physical domain (z,r) to the computational 
domain (^,T|) where 


% = 4(z.r) 

and 

r\ = Ti(z,r) 

Applying the chain rule of partial differentiation to this pair yields 

B_ = B^B_ 

dz dz 0 ^ Bz 

b _ = as a_ an _d_ 

dr dr 3^ dr dn 
and 

d£= ^dz + ^ dr 
dn = nz dz + nr dr 


So 


Similarly, 


d£ 



[dzl 

.dn. 


.Vzllr. 

LdrJ 


dz 


'*5^1 

di 

dr. 


j 1 

1 

dn. 
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Therefore, 




’ Z $V 

_Tlz Tlr. 


A r n. 


Taking the inverse yields 


Tlz Tlr 



- z n 


where J is the Jacobian of the transformation given by 


J — 4z ^Ir " Vz 

= l/(z^ rT1 - Zn r^) 

= 1 / Vy = 1 / (cell volume) 


and the metrics are 


\z - J r Tl 

rf 

Hi 

1 

II 

Tlz = - J ^ 

“Hr = J ^ 


By defining the normalized metrics 


kz 

It 1 

^ 


v^ 2 + ^ 2 

Vzt , 2 + r^ 2 

■ V*/ + s, 2 

Vz ,! 2 + r-q 2 

% 

* r $ 

, it - 


VtJz 2 + T| r 2 

Vz ^ 2 + r ^ 2 

VrU 2 + T| r 2 

Vz ^ 2 + r ^ 2 


the normal and tangential velocity components (see Fig. B.l) to surfaces of constant £, and 
T| can be defined [46 and/or 47]. 
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T| = CSt 


V* 



■ 

i 

i 

£ = cst 


Figure B.l Normal and Tangential Velocity Components at 
Surfaces of Constant tj and ^ 

The velocity components normal and tangential to a surface of ^ = cst. located on 
the surface at (i + 1/2, j) are 

u’ = sj z u + s^ v 
_ u + (-zn) v 
d A 

where d A = V r^ 2 + z^ 2 

The velocity components normal and tangential to a surface of t) = cst. located on the 
surface at (i, j + 1/2) are 

V' = Sjz U + Sjr V U' = Sjr U - Sj z V 

-r^u + z^v z^u - (-r^)v 

da de 

where dfi = V r^ 2 + z^ 2 

Alternatively, the contravariant velocity components can be defined as 



U = + £rV 

V = T| Z U + TlrV 

Therefore, 

LL = u - dA 

f “ VdB 


v’ = -Sjr u + Si z V 

= + r n v 

ti A 
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Applying this transformation to the governing equations 


dQ dF 1 
— + 

dt dz 


BG' 
B r 


+ H' = 0 


yields 

Qt + + tIzFt, + + T| r G 7 ^ + ^=0 

Dividing by the Jacobian and adding and subtracting like terms [27] yields 


lQ\ k ^rG' \ (tlzFjjT M £ 

W J k ( J h 3 



(Bl) 

(B2) 


(B3) 


Substituting the metrics into the last two terms of Eq. (B3) and performing the 
indicated differentiation (switching the order of differentiation first) shows these terms sum 
to zero. Defining the quantities 




F = 


£z F' + £ r G' 


G = 


rizF + Tj r G' 


(B4) 



and using the definition of the contravariant velocity components, the governing equations 
in the transformed coordinate system may be written 


dQ 

dt 


+ 


dF ’ dG ' 

+ + 

d£ dq 


H’ = 0 


(B5) 


where the vectors Q, F', G’, H* are the Q, F, G\ H' vectors in the transformed 
coordinate system. Splitting the F, G\ and H' vectors into inviscid and viscous 
components yields 


dQ dF df; 

— + — + 

dt B% d$ 


+ 



dn 


dG v a — 
— + H + Hv= 0 
dq 


(B6) 
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Next the Thin Layer Approximation (TLA) is applied to the governing equations. This 
assumes viscous terms with streamwise derivatives are small compared to the viscous 
terms containing transverse derivatives, and causes terms containing 


to be neglected. Thus, the term 


d_ 

dEy 


is neglected and the partial derivatives become 

a _ a 


dz 112 an 


and 


dr dr] 


The shear stress and heat conduction terms become 

.1 


trr = 

3 


_ dv du 

2rir Tlz — 

dri dq 


T ee - 


= -2i 



dv 
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n r — 

+ nz— 


an 

an l 


* 


- in y. 


n f 


Cl = - 2 -n 
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_ du dv 
2 T]z ~ tlr 

dri dr] 


^ * 


~ V- 




trz = - M- 


n r - 


Tlz- 


9w 

dn. 

dw 

dnj 
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+ H ? 


du dv 

^dn +Tlz dn 
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dT 9T t 9T dT 

— , — => \z — . "Hr — 

3z 9r 9t) dr| 


Using the rotated velocity components u’ and v' the final set of governing equations 
and components are 


L*± + *l 

J at a^ 


3G aG v 

+ — + — 

an arj 


+ H + H v = 0 


where 


(B8) 



F = d A 


pu 

puu' + 4-r^P 
“A 

pvu' + ^-z^P 
dA 
pwu' 

(El + p) u’ 


G 


d B 


pv’ 

puv' + 

dB 

pvv' + 

d B 

pwv* 
(E, + p) 


v’ 


where u' and v' are the rotated velocity components normal to surfaces of constant £ and T|, 
respectively. 
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0 


Hv = 


- » - rkk u ’i- r « v i) 

+ rk (j? - 


rs?(j7 -^ Wi >) 


k 


Rl_ 


rRePr(y 0 -l) J Re P t(y<> - 1) <* 


J Re Pr (Yo - 1) Sr 2RePr(Yo- l)^ jr * T ^ Ty] 


It should be noted that j r , j z , B e , represent the radial current component, axial current 
component, and circumferential magnetic field, respectively and are not the radial, axial, 
and circumferential derivatives. All other subscripts denote a derivative with respect to the 
subscripted variable. 

Lastly, using the chain rule of partial differentiation, the coordinate transformation 
is applied to the equation governing the electromagnetic field to yield 
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Appendix C 


Derivation of the 
Inviscid Flux Vector Jacobians 


Flux splitting of the rotated inviscid flux vectors, F and G, is achieved by 
decomposing the inviscid flux into two components. One component contains the 
eigenvalues traveling in the positive coordinate direction and the other component contains 
the eigenvalues traveling in the negative coordinate direction. The purpose of this appendix 
is to illustrate how the Jacobian matrices of F and G are diagonalized by the matrix of 
eigenvalues, X. 

The rotated inviscid flux vectors can be represented as the product of a rotated 
Jacobian matrix, A and B, respectively, and the solution vector Q 

F = f^Q=AQ G = ^Q=BQ (Cl) 

dQ dQ 

Recalling the definition of F from Appendix B, for example 


and 


F = d A 



pu' 


puu' 

+ 

£ r " p 

pvu' 

+ 



pwu' 


l E ' 

+ p) 

u’ 


u' = SizU + SirV 


where u' is the velocity normal to a surface of £= cst and sj z and Sjj- are the normalized 
metrics presented in Appendix B. The solution vector is given by 


a 


144 




P 


roti 


pu 


02 

Q = 

pv 

= 

03 


pw 


04 


. Et . 


LosJ 


The rotated Jacobian matrices can also be written in terms of the vector of primitive 
variables, V, 


and 


3F 

_9Q 

aV dF 

av 

aq 

av 

dQ av 

50 

dG 

= ao. 

avaS 

av 

5Q 

av 

dQ av 

dQ 


(Cl) 


where 


V = [ p u v w P ] T = [Vi V 2 V 3 V 4 V 5 ] t 


Define 


3Q 


and 


s-i-32 

av 


(C2) 


0V 

and diagonalize the Jacobian — (as opposed to — , which is more difficult to 

3Q av dQ 

diagonalize) as 

cx>a aCa 

dQ av 


The rotation matrix, R A , which contains the normalized metrics and which transforms the 
velocities between the Cartesian and the transformed coordinate system, can be defined 
[refer to 29, 60] so that C A - = C A Ra and the Jacobian becomes 


av 3F 
3Q av 


= R a 1 C a 1 A a C a R a 


(C3) 
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where 


Ra 


1 0 0 0 0 

0 Sji Sjj 0 0 

0 ■ Sjj Sij 0 0 

0 0 0 1 0 

.0 0 0 0 1 - 


Substituting Eqs. (C2) and (C3) into (Cl) yields 

A = S'^C^Aa C a R A S 

Similarly (C4) 

B = S* 1 Rb 1 Cb 1 A b C b Rb S 


The Jacobian S is found by first writing the vector of primitive variables, V, in 
terms of the solution vector, Q. 


Qi 

02 
Qi 
Ql 
Qi 

Qt 
Qi 

From which the following may be evaluated 
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V 4 
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8Vi 8Vi 9Vi 9Vi 9Vi 
9Qi dQ2 ^03 aQ 4 3 Qs 
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av s av 5 

dQi dQs _ 
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Thus, 


S = 


1 0 0 0 0 

-u/p 1/p 0 0 0 

-v/p 0 1/p 0 0 

-w/p 0 0 1/p 0 

P a -P u *P v -P w P 


a = |-(u 2 + v 2 + w 2 ] 


where P = (y - 1 ) ”2 

In a similar manner S' 1 is found, by placing Q in terms of V and evaluating 
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0 
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0 

0 0 
P 0 


w 0 

a p u p v p w 1/P 

Examination of Eq. (C3) shows one unknown matrix on the left hand side of the 

dF ^ — • 

expression, — , which can be evaluated, by writing F in terms of V, to be 

8V 
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au'd A 
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-u'd A 
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Rearranging Eq. (C3) and using the definition of S we can write 


C A J A a C a = 


R A S 



(C5) 


All matrices on the right hand side of Eq. (C5) are known. Using MACSYMA (an 
interactive symbolic algebra program from Symbolics, Inc., Eleven Cambridge Center, 
Cambridge, MA 02142) the eigenvalue problem posed by Eq. (C5) was solved. C A * is 
composed of columns of right eigenvectors and A a is a diagonal matrix of eigenvalues. C A , 
as opposed to C A \ is given by 


where 


p-i 

'■'A 


10 0 11 
0 0 0 c/p -c/p 

0 10 0 0 

0 0 10 0 

0 0 0 c 2 c 2 


Aa 


d A u’ 1 

d A u’ 0 

d A u’ 

0 d A (u’ + c) 

d A (u’ - c) _ 


c 2 = 


11 

P 


The matrices diagonalizing B, Rb, Cb. and Ab, can be obtained directly by 
substituting 


v‘ for u’ 


and 


in R a , C a , and A a . 


Sir* 


da for d A 
-Sj z for Sj z , Sir 


148 



Appendix D 


Implicit Viscous and Source Term Matrices 


The matrices M a and M5 arise from the derivation of the implicit differencing of the 

/s 

viscous flux vector, 5 G v ,where 


G v = 


Re 
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I 

J (P4 Utj + p 2 V^j ~ ^^V 

Oi 1 

J (P2 Ur! + pl Vfjj + V 

J p3 ^ Z5 W 

1 b,9t 

L Re Pr(y 0 - lj ’ 


By partitioning the vector into: 1 ) the terms common to both the cartesian and 
cylindrical formulations and 2) the additional terms arising from the cylindrical 

/s 

axisymmetric formulation ( in G v shown above, the three terms containing 1/r are the 
additional terms arising from the cylindrical coordinate formulation) such that 

/s /\ 

Gy = Gvcom + Gvcyl 

Differentiating this expression with respect to time, multiplying by At and using the 
/\ 

definition of 8Fy 


so, ■ or' - a; = 

ot 


we can write 


5 G V = 8G 

vcom + 5 G 


*vcyl 
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Employing the vector of primitive variables, V 

V = (p u v w T) t 


we can write 


5G V = M.|-5V + 5G V cy i 

dr] 


because 5G V contains derivatives of the primitive variables with respect to T| only. The 

*\ 0 -* xv 

quantity 5V can be derived and then M a found such that M a -^5V = 8G V . M a is 
found to be 


M a = 
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•**(3 A * ’ll 
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p- - •(‘H + f *1) 


p 5 = -k(4 + >|) 


We can change variables from 8V to 8Q by defining the Jacobian N where 
jsj _ d Y. suc h that 8V = N 8Q. In addition, 8G V cyi can be expressed in terms of 8Q such 


3Q 


that 


8G V = + M b 8Q 
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The Jacobian, N, is found to be 


where 
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N =J- 
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- v 

- w 


|_e(2a - e t ) 
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- e u 


0 0 0 
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by evaluating 


avi avi avi avi avi 
dQi 9Q2 dQ3 ^Q4 dQs 
I - - - I 


N = 


av 5 av 5 

_9Qi “ “ “ 305 


after placing V 1-V5 in terms of Q1-Q5. It should also be noted that in deriving the fifth row 
of N, y was assumed to be locally independent of Q. 

In a similar manner G v C yi is expressed in terms of Q 1 -Q 5 and then M5 is obtained 
by evaluating 
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The implicit in viscid : 

source term, 8Hj 

ij, is placed in 

terms of 8Q by writing 


8Hi,j = Qij 5Qij 


151 



The inviscid source term matrix, Cj is found by rewriting H, 



pv 
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puv 

r 



pw 2 

r 



in terms of QrQs, and then evaluating Cj = — — to obtain 

3Q 


where 
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a _ u 2 + v 2 4- w 2 
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It should also be noted that in deriving the fifth row of Q, as was done for N, y was „ 
assumed to be locally independent of Q. 
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